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Abstract 

The  wealth  of  available  scientific  data  collected  by  the  modem  constellation  of 
meteorological  satellites  can  be  exploited  in  new  and  innovative  ways,  with  direct  benefit 
to  directed  energy  applications.  The  Atmospheric  Infrared  Sounder  (AIRS)  provides 
accurate,  geolocated  temperature  data  for  100  vertical  levels,  with  near  global  coverage 
twice  daily.  AIRS  temperature  data  are  used  as  the  starting  point  for  the  derivation  of 
atmospheric  parameters  critical  to  prediction  of  sensor  or  directed  energy  system 
performance.  The  accuracy  and  utility  of  the  atmospheric  characterization  can  be  improved 
by  combining  the  data  from  multiple  sources.  Twenty-four  hour  data  coverage  is  achieved 
by  extending  the  atmospheric  characterization  between  satellite  observations  using 
Numerical  Weather  Prediction  (NWP)  models.  A  technique  is  developed  to  derive  wind 
profiles  using  AIRS  temperature  data  and  the  accuracy  of  these  winds  is  evaluated.  These 
winds  are  not  found  to  be  as  accurate  as  the  winds  available  from  modern  NWP  models, 
but  their  usefulness  is  demonstrated  for  the  stratosphere  and  lower  mesosphere  where  NWP 
model  data  are  not  readily  available.  Techniques  to  derive  optical  turbulence  strength,  as 
modeled  by  index  of  the  refraction  structure  function  C2,  are  applied  to  integrated 
meteorological  satellite  data  and  NWP  data.  The  existing  technique,  which  bases  the 
derivation  of  Cn2  for  optical  wavelengths  on  vertical  atmospheric  temperature  gradients,  is 
improved  by  including  additional  wind  shear  components  and  local  pressure  gradient 
deviations  from  the  traditionally  assumed  hydrostatic  balance. 
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OPERATIONAL  EXPLOITATION  OF  SATELLITE-BASED  SOUNDING  DATA 
AND  NUMERICAL  WEATHER  PREDICTION  MODELS  FOR  DIRECTED 

ENERGY  APPLICATIONS 

1  Introduction 

The  performance  of  directed  energy  and  imaging  systems  is  affected  by  the  state  of 
the  atmosphere  over  the  path  of  interest.  Temperature  gradients,  optical  turbulence,  and 
the  presence  of  aerosols  and  clouds  all  have  significant  impact,  but  these  quantities  can  be 
difficult  to  accurately  measure  or  predict  over  extended  paths  through  the  atmosphere.  This 
research  examines  the  availability,  resolution,  and  accuracy  of  data  collected  by  several 
satellite-based  sounding  instruments,  and  analyzes  their  applicability  to  new  atmospheric 
characterization  methods  that  will  benefit  the  directed  energy  community.  The  advantages 
of  combining  satellite-based  sounding  data  with  the  output  from  modern  Numerical 
Weather  Prediction  (NWP)  models  is  detailed  and  innovative  techniques  to  improve  the 
methods  used  to  derive  the  strength  of  atmospheric  optical  turbulence  are  developed. 

1.1  Motivation  for  this  Research 

Optical  turbulence  distorts  and  spreads  electromagnetic  (primarily  in  the  visible  and 
near-IR)  radiation  propagating  through  the  atmosphere,  and  the  ability  to  characterize  the 
strength  of  the  optical  turbulence  along  a  specified  slant-range  path  using  data  from  recent 
satellite  measurements  would  have  both  civilian  and  military  operational  value.  In  order 
to  accurately  predict  the  effect  that  the  atmosphere  will  have  on  a  sensor  or  directed  energy 
system,  it  is  necessary  to  be  able  to  characterize  the  path-weighted  strength  of  the  optical 
turbulence  and  describe  the  properties  of  any  clouds  or  aerosols  along  the  path.  When 
designing  directed  energy  systems  and  predicting  atmospheric  effects,  the  common 
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approach  is  to  make  assumptions  about  the  state  of  the  atmosphere,  applying  a  simple 
atmospheric  model  or  long-term  average  values.  This  can  lead  to  unrealistic  expectations 
about  how  a  particular  directed  energy  system  would  perform,  or  a  lack  of  appreciation  for 
the  variability  of  performance  as  atmospheric  conditions  change. 

Atmospheric  turbulence  can  be  measured  using  a  thermosonde,  but  this  method  has 
limited  usefulness  for  directed  energy  applications.  The  thermosonde  measures  a  profile 
above  a  single  geographic  location,  which  does  not  provide  the  information  needed  to 
characterize  an  extended  slant-range  path  from  platform  to  target.  A  scintillometer  can 
measure  optical  turbulence  along  a  path,  but  it  also  has  limited  utility  for  long-range 
directed  energy  applications.  Two  large  drawbacks  to  operational  use  of  a  scintillometer 
for  directed  energy  applications  are  the  requirement  to  have  equipment  at  both  ends  of  the 
path,  and  the  possibility  that  for  a  long  path  with  strong  turbulence,  saturation  of  the 
scintillometer  makes  accurate  measurement  impossible. 

There  is  a  wealth  of  space-based  atmospheric  data  that  can  be  used  to  derive  the 
required  atmospheric  characterization.  The  continuous  stream  of  data  collection  by  dozens 
of  scientific  instruments  provides  the  basis  for  both  improvements  in  our  atmospheric 
models  and  derivation  of  specific  parameters  of  interest.  Starting  with  satellite-measured 
temperature  data,  it  is  possible  to  derive  3-D  fields  of  winds,  gradient  Richardson  numbers, 
and  optical  turbulence  strengths,  as  characterized  by  the  index  of  refraction  structure 
function,  C,,2. 

The  focus  of  this  research  is  to  analyze  data  collected  by  various  scientific 
instruments  in  operation  on  meteorological  satellites  and  to  develop  and  evaluate 
techniques  to  take  full  advantage  of  the  data  these  instruments  can  provide,  and  show 
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innovative  ways  the  data  from  different  sources  can  be  combined  to  achieve  operationally 
useful  atmospheric  characterization.  The  types  of  space-based  instruments  studied  include 
IR  and  microwave  sounders,  visible/IR  imagers,  and  LIDAR  (Light  Detection  and 
Ranging)  units.  The  primary  objective  is  to  demonstrate  the  ability  to  create  operationally 
relevant,  near  real-time,  3-D  gridded  data  sets  of  atmospheric  turbulence  and  cloud 
characterizations,  derived  from  satellite  observations  and  NWP  models. 

1.2  Document  Structure 

This  document  is  organized  into  five  main  chapters.  This  first  chapter  provides  an 
introduction  to  the  research  goals  and  the  motivation  for  the  project,  along  with  a  general 
description  of  the  operational  capabilities  that  this  effort  enables.  The  second  chapter 
details  the  background  scientific  basis  for  the  techniques  employed  and  includes  a  survey 
of  the  satellite  instruments  and  data  sources  investigated.  The  relative  strengths  and 
weaknesses  of  the  data  provided  by  the  individual  instruments  and  weather  models  are 
discussed,  as  is  the  role  each  data  source  can  play  in  the  combined  atmospheric 
characterization  techniques.  This  discussion  is  followed  by  a  review  of  the  topic  of 
atmospheric  turbulence,  covering  the  background  information  necessary  for  the 
development  of  modified  turbulence  calculation  techniques  for  a  well-mixed  atmospheric 
layer. 

The  third  chapter  covers  the  methodology  for  the  research.  Several  new  or 
expanded  applications  of  satellite-derived  applications  are  examined  in  this  work.  The 
three  primary  tasks  were  analyzing  the  height  dependent  error  and  bias  present  in  the 
satellite-sounding  data  used  in  this  research,  expanding,  improving,  and  evaluating  the 
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satellite-based  windfinding  algorithm  originally  developed  in  the  authors  2010  master’s 
thesis  [1],  and  developing  several  techniques  to  derive  optical  turbulence  strength  from  the 
sounding  and  NWP  data.  These  three  topics  are  detailed  in  Sections  3.1,  3.2,  and  3.3.  The 
methodology  for  the  additional  topics  of  satellite-based  cloud  location,  cloud-top  height 
determination,  and  techniques  to  use  the  NWP  model  data  to  fill  in  temporal  gaps  between 
satellite  overpasses  is  also  presented  in  Chapter  3. 

Chapter  4  presents  and  analyzes  the  results  obtained  by  implementing  the 
techniques  developed  and  explained  in  Chapter  3.  These  two  chapters  share  a  parallel 
structure,  with  Sections  4.1,  4.2,  and  4.3  detailing  the  results  from  the  methods  described 
in  Sections  3.1,  3.2,  and  3.3.  Chapter  5  summarizes  the  impact  of  the  project,  details  the 
conclusions  drawn  from  this  research  endeavor,  and  outlines  recommendations  for  areas  of 
follow-on  research. 

The  satellite  sounder  and  imager  data  used  throughout  this  research  are  collected 
by  instruments  onboard  satellites  that  are  part  of  NASA’s  Earth  Observing  System  (EOS). 
These  data  are  transmitted  to  ground  stations  and  relayed  to  White  Sands,  NM  via  the 
Tracking  and  Data  Relay  Satellite  System  (TDRSS).  The  data  are  then  sent  to  the 
appropriate  Distributed  Active  Archive  Center  (DAAC)  for  additional  processing  and 
distribution.  All  of  the  satellite  data  sets  used  in  this  research  are  publically  available  for 
download  from  either  the  Goddard  Earth  Sciences  Data  and  Information  Service  Center 
(GES-DISC)  or  the  EOS  Data  and  Operations  System  (EDOS). 
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2  Background  and  Literature  Review 


2.1  Infrared  Sounder  Theory 


2.1.1  Blackbody  Theory 

The  wavelength  dependent  intensity  of  the  upwelling  radiation  detected  by  a 
meteorological  satellite  is  dependent  on  the  temperature  of  the  atmosphere  as  described  by 
Planck’s  function  [2] 


Bxcn  = 


2  he 


A,5  (ehc,k°AT -l) 


(1) 


where  c  =  2.998  x  108  m/s  is  the  speed  of  light,  h  =  6.626  x  10"34  J  s  is  Planck’s  constant 
and  Jcb  =  1.381  x  10'23  J/K  is  Boltzmann’s  constant. 

The  wavelength  at  which  the  Planck  function  peaks  indicates  the  temperature  of  the 
blackbody,  according  to  Wien’s  Displacement  Law  [3] 

4,„  -  ^  (2) 


where  the  constant  kw  =  2897. 9//m  K  . 

To  use  these  relationships,  it  is  necessary  to  determine  how  closely  the  atmosphere 
approximates  a  perfect  blackbody  in  the  infrared.  The  emittance  of  a  body  can  be 
wavelength  dependent  and  is  defined  as: 

emitted  •  radiation  -at- A, 


BAT) 


(3) 


A  perfect  blackbody  has  an  emittance  equal  to  one,  and  by  Kirchhoff  s  Law,  would  have 
an  absorptance  equal  to  one  if  it  is  in  local  thermodynamic  equilibrium.  Below  about  100 
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km  altitude  in  Earth’s  atmosphere,  collisions  between  molecules  are  frequent  enough  to 
make  the  assumption  of  local  thermodynamic  equilibrium  reasonable  [3].  For  satellite 
meteorology  applications,  it  is  a  good  approximation  to  say  that  wavelength  dependent 
absorptance  equals  the  wavelength  dependent  emission,  ax=ex. 

2.1.2  Satellite  Retrieval  of  Temperature 

The  ability  to  look  down  from  space  and  measure  atmospheric  temperatures  at 
multiple  depths  within  the  atmosphere  relies  on  molecular  absorption  by  atmospheric 
gases.  By  knowing  the  concentration  of  a  gas  in  the  atmosphere  and  its  mass  absorption 
coefficient  for  a  given  wavelength,  an  optical  depth  can  be  determined.  When  a  space- 
based  sensor  with  a  narrow  field-of-view  collects  upwelling  radiation  within  a  narrow 
wavelength  channel,  the  altitude  from  which  most  of  that  radiation  emanated  can  be 
determined. 

If  one  of  the  constituent  molecules  that  make  up  the  atmosphere  strongly  absorbs 
the  wavelength  being  measured,  any  emissions  from  the  surface  or  lower  atmosphere  will 
not  make  it  to  space.  Instead,  almost  all  of  the  emission  received  will  be  from  an  altitude 
near  the  top  of  the  atmosphere.  For  a  different  channel,  one  with  a  wavelength  for  which 
the  atmosphere  is  transparent,  most  of  the  radiance  measured  by  a  space-based  sensor 
emanates  from  the  warmest,  thickest  layer  of  the  atmosphere  (the  very  bottom  of  the 
atmosphere)  and  the  surface  of  the  Earth.  Between  these  two  extremes,  each  channel  falls 
somewhere  along  the  sloped  edge  of  a  molecular  absorption  line,  and  a  weighting  function 
is  calculated,  as  described  in  Equation  (4)  and  Figure  1 .  The  height  at  which  the  weighting 
function  peaks  is  important.  The  sensor  receives  a  larger  fraction  of  the  atmospheric 
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thermal  emission  from  that  height  than  from  any  height  above  or  below  it.  Schwarzchild’s 
equation  is  manipulated  to  yield  a  weighting  function  for  each  individual  wavelength. 

d  f  V\ 

wAh’M)  =  —  (4) 

dh  V 

y 

where  Ms  a  generalized  height  coordinate,  ju  is  the  cosine  of  the  zenith  angle,  and  z  }y‘ 

is  the  slant  path  transmittance  from  level  h  to  the  satellite.  For  each  wavelength  measured, 
the  peak  in  this  weighting  function  is  used  to  assign  a  height  to  the  radiance  measured.  The 
process  of  assigning  temperatures  to  these  heights  is  more  complicated — an  iterative 
process  of  solving  to  find  a  full  temperature  profile  that  would  result  in  the  radiance 
measured  across  all  of  the  sounder’s  channels  [3]  [4], 

The  most  commonly  used  gases  in  temperature  sounding  are  carbon  dioxide,  water 
vapor,  oxygen  and  ozone  (carbon  dioxide  and  oxygen  having  the  advantage  that  they  are 
well-mixed  to  an  altitude  of  80  km,  making  it  possible  to  determine  their  densities  with 
altitude)  [5].  Without  knowing  concentration,  the  altitude  for  water  vapor  returns  is 
difficult  to  determine,  but  water  vapor  channels  are  useful  for  imagery.  Common 
absorption  wavelengths  used  by  satellites  are  in  the  vicinity  of  15  pm  (carbon  dioxide),  9.6 
pm  (ozone),  and  5-8  pm  (water  vapor)  [3].  The  primary  satellite  instrument  exploited  in 
the  research,  the  Atmospheric  Infrared  Sounder  (AIRS),  additionally  takes  advantage  of 
the  15  pm  and  4.3  pm  carbon  dioxide  absorption  features  [6]. 

The  determination  of  temperature  by  passive  atmospheric  sounding  starts  with 
Schwarzchild’s  equation,  shown  in  Equation  (5).  Vertical  sounding  theory  relies  on  the 
integrated  form  of  this  equation: 
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4=V"'+JV  "  Bt(T) 


<5q-<5a 


dS± 

A 


(5) 


Lx  is  the  radiance  measured  by  the  satellite  instrument,  is  the  radiance  incident  on  the 


absorber  of  interest,  Bx  ( T )  is  the  Planck  function,  and  //  is  the  cosine  of  the  zenith  angle. 


S0  is  defined  as  the  vertical  optical  depth  of  the  entire  atmosphere  and  the  wavelength 
dependent  optical  depth  Sx  is  used  to  determine  height  in  the  atmosphere. 

The  usefulness  of  an  individual  measurement  is  limited,  but  when  measurements 
are  combined  from  a  series  of  closely  spaced  channels,  each  receiving  a  slightly  different 
wavelength  near  a  strong  molecular  absorption  line,  a  temperature  profile  can  be 
constructed.  Since  the  gas  has  a  different  optical  depth  for  each  wavelength,  each  channel 
of  the  sensor  will  be  most  sensitive  to  radiance  emitted  by  a  different  altitude  in  the 
atmosphere.  Matching  these  individual  measurements  to  a  first-guess  model  of  the 
atmosphere  through  an  iterative  process  can  yield  a  full  temperature  profile. 

The  radius  of  typical  operational  polar  orbits  are  much  lower  than  geostationary 
orbits  (700-800  km  vs.  35,786  km),  and  therefore  do  not  need  short  wavelengths  to  have 
useful  resolution  [3].  For  this  reason,  instruments  on  a  polar,  low  Earth  orbiting  satellite 
can  take  advantage  of  microwave  absorption  features  like  oxygen’s  strong  rotation  feature 
at  60  GHz  (5  mm).  The  Advanced  Microwave  Sounding  Unit  (AMSU)  is  an  instrument 
employed  on  polar  orbiting  satellites  which  uses  15-20  channels,  several  on  the  edge  of 
this  oxygen  absorption  band. 

The  weighting  functions  for  these  AMSU  channels  plotted  against  pressure  in 
millibars  (mb)  is  shown  in  Figure  1.  One  millibar  is  equal  to  one  hectopascal  (hPa).  For 
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each  wavelength,  approximate  altitude  can  be  determined  using  the  pressure  for  the  point 
where  the  corresponding  channel’s  weighting  function  reaches  a  maximum  [2], 


0  .1  2  3  A  .5  .6  .7 
WEIGHTING  FUNCTIONS 


Figure  1.  Microwave  weighting  functions  for  AMSU-A  sounding  channels 
(Courtesy  NASA/GES-DISC)  [7]. 

2.2  The  Atmospheric  Infrared  Sounder  (AIRS) 

High  resolution  infrared  sounders,  such  as  the  AIRS  instrument  onboard  NASA’s 
Aqua  satellite,  provide  global  atmospheric  temperature  measurements  with  excellent 
vertical  resolution,  comparable  to  the  measurements  made  by  weather  balloons.  AIRS  is 
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an  infrared  spectrometer  with  2378  channels  covering  a  spectral  range  between  3.7  and 
15.4  pm,  taking  advantage  of  carbon  dioxide  absorption  features  at  both  4.3  pm  and  15  pm. 
The  Aqua  satellite’s  polar  orbit  and  the  AIRS  instrument’s  swath  width  of  about  1600  km, 
lead  to  soundings  covering  most  of  the  planet  twice  each  day.  The  IR  instrument  does  not 
retrieve  measurements  from  below  the  highest  layer  of  clouds,  but  microwave  soundings 
from  the  Advanced  Microwave  Sounding  Unit  (AMSU),  also  onboard  Aqua,  are  used  to 
fill  in  missing  surface  and  lower  atmosphere  temperature  measurements  [1].  Figure  2 
depicts  the  scan  geometry  for  the  AIRS  instrument,  detailing  field-of-view  geometry  and 
scan  pattern  overlap  [8]. 
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•  Altitude:  705  km 
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•  Ground  Footprints:  90/Scan 


Figure  2.  Scan  geometry  for  the  AIRS  instrument  (Courtesy  NASA/GES  DISC)  [9]. 

During  each  99  minute  orbit  of  the  satellite,  the  Earth  rotates  approximately  25 

degrees  about  its  axis,  so  that  each  Aqua  ground  track  crosses  the  equator  approximately 

1730  km  west  of  the  previous  crossing,  as  shown  in  Figure  3.  Near  the  equator,  the 

coverage  is  not  complete  and  the  surface  area  missed  is  filled  in  by  the  next  day’s  flight 
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path.  At  latitudes  greater  than  about  45  degrees,  each  swath  overlaps  the  previous  one 
causing  some  locations  to  be  sounded  four  (or  more)  times  per  day.  Most  locations  on 
Earth  fall  within  the  AIRS  data  swath  twice  in  a  24  hour  period,  under  both  an  ascending 
pass  (during  the  local  afternoon)  and  a  descending  satellite  pass  (during  the  local  early 
morning).  One  of  the  primary  objectives  of  this  research  is  the  development  and  evaluation 
of  new  techniques  to  exploit  this  AIRS  temperature  data  for  directed  energy  applications. 


Figure  3.  AIRS  coverage  during  three  sequential  ascending  passes  on  1  Oct  2015. 
The  left  image  is  the  first  pass  between  1726  and  1816  UTC.  The  center  image 
shows  the  first  pass  and  the  second  pass,  between  1906  and  1956  UTC.  The  right 
image  shows  the  first  two  passes,  and  adds  the  third  2046  to  2136  UTC  pass. 

2.3  Additional  Meteorological  Satellite  Instruments  Investigated 

The  AIRS  instrument  has  excellent  vertical  resolution,  but  poor  horizontal  resolution 
when  compared  to  imaging  instruments,  and  it  has  the  limitation  of  not  being  able  to 
retrieve  temperatures  below  thick  cloud  cover.  AIRS  is  one  of  five  scientific  instruments 
on  the  Aqua  satellite,  and  Aqua  is  part  of  NASA’s  A-Train  constellation  of  Earth-observing 
satellites.  These  satellites  follow  the  same  orbital  track  with  the  trail  satellite  only  a  couple 
minutes  behind  the  lead.  The  common  ground  track  and  near  simultaneous  observations 
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of  each  satellite  instrument  enable  multiple  data  sets  to  be  used  together,  so  that  the 
strengths  of  one  type  of  instrument  can  make  up  for  the  shortcomings  of  another. 

The  orbits  of  the  satellites  in  the  A-Train  constellation  are  sun-synchronous, 
meaning  the  orbit’s  98.14°  inclination  is  chosen  so  that  for  its  semi-major  axis  of  7078  km, 
the  right  ascension  of  the  orbit’s  ascending  node  precesses  at  the  same  rate  as  the  Earth’s 
revolution  around  the  sun.  This  causes  the  ground  track  to  pass  over  each  location  on  the 
Earth  at  approximately  the  same  local  time  every  day,  throughout  the  year  [3].  The  A- 
Train  gets  its  name  from  the  fact  that  the  satellite  formation  always  crosses  the  equator 
northbound  in  the  early  afternoon  local  time  for  the  region  directly  below — the  ‘Afternoon 
Train.’  Figure  4  shows  the  satellites  and  instruments  that  make  up  the  NASA  A-Train. 


Figure  4.  NASA  depiction  of  the  A-Train  constellation  showing  sensor  geometry  for 
the  primary  Earth-observing  instruments  (Courtesy  NASA/JPL-Caltech)  [10]. 
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2.3.1  Advanced  Microwave  Sounder  (AMSU) 

The  AMSU  began  operation  in  1998,  and  is  onboard  satellites  NOAA-15  through 
NOAA-19,  NASA  Aqua,  and  EUMETSAT  MetOp  A  and  MetOp  B.  The  instrument 
consists  of  two  separate  sounders.  These  sounders  are  the  AMSU-A  (with  15  channels 
between  23.8  and  89  GHz)  and  the  AMSU-B  (with  5  channels  between  89.9  and  183.3 
GHz).  The  AMSU  employs  a  whisk  broom  scan  pattern,  sweeping  49.5  degrees  on  either 
side  of  nadir  for  a  1650  km  wide  swath.  Nadir  is  the  viewing  angle  of  the  centerline  of  the 
data  swath,  directly  below  the  satellite’s  track.  During  each  6  second  sweep,  radiance 
values  for  every  channel  are  collected  for  30  fields  of  view  for  the  AMSU-A.  The  higher 
frequency  AMSU-B  collects  radiance  from  90  fields  during  each  sweep.  The  cross-track 
resolution  at  nadir  is  45  km  for  AMSU-A  and  15  km  for  AMSU-B  [3]  [8]. 

The  AMSU  sounders  work  well  in  combination  with  infrared  sounders,  which  are 
not  useful  for  tropospheric  measurements  below  cloud  cover.  Because  the  AMSU  has  only 
20  channels,  the  vertical  resolution  of  the  sounding  data  does  not  match  the  infrared 
sounders,  but  the  microwave  frequencies  are  able  to  pass  through  cloud  cover  and  measure 
brightness  temperatures  all  the  way  to  the  surface.  The  AMSU  senses  upwelling  terrestrial 
radiation,  comprised  of  a  combination  of  radiation  emitted  by  Earth’s  surface  and  radiation 
emitted  by  the  atmosphere.  As  discussed  in  Section  2.1.2,  each  channel  is  radiatively 
selective,  primarily  detecting  microwave  radiation  from  a  discrete  layer  within  the 
atmosphere.  The  height  of  this  layer  corresponds  to  the  height  at  which  that  channel’s 
weighting  function  reaches  a  maximum  value.  On  NASA’s  Aqua  satellite,  the  standard 
AMSU-B  was  replaced  with  a  modified,  4-channel  version,  called  the  Humidity  Sounder 
for  Brazil  (HSB).  The  HSB  stopped  operating  in  February  of  2003  due  to  an  electrical 
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failure  within  the  components  that  control  the  scanning  mirror,  thus  the  only  microwave 
sounding  data  available  from  the  Aqua  satellite  since  that  time  are  for  the  AMSU-A 
channels. 

2.3.2  Moderate  Resolution  Imaging  Spectroradiometer  (MODIS) 

The  Moderate-resolution  Imaging  Spectroradiometer  is  a  scientific  instrument 
onboard  NASA’s  Earth  observing  Aqua  and  Terra  satellites.  MODIS  captures  radiation  in 
36  spectral  bands  with  wavelengths  between  0.4  and  14.4  pm.  The  spacing  of  these 
spectral  bands  relative  to  atmospheric  absorption  features  is  shown  in  Figure  5.  An 
advantage  of  the  MODIS  data  is  the  spatial  resolution  for  these  channels,  providing 
horizontal  resolution  ranging  from  250  m  to  1  km. 


2100 


2500  2700 


WAVELENGTH  (nm) 


Figure  5.  MODIS  channel  wavelengths  (Courtesy  NASA)  [11]. 
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In-flight  calibration  of  the  data  is  accomplished  using  an  internal  blackbody,  a  solar 
diffuser,  and  a  spectral  radiometric  calibration  assembly.  Each  MODIS  instrument  images 
a  2330  km  wide  swath,  collecting  data  for  almost  all  of  the  Earth’s  surface  twice  daily. 
Figure  6  shows  the  coverage  of  24  hours  of  sequential  descending  MODIS  passes. 


Figure  6.  MODIS  imagery  from  the  14  descending  passes  during  one  24  hour 
period.  This  image  was  retrieved  from  https://lpdaac.usgs.gov,  maintained  by  the 
NASA  EOSDIS  Land  Processes  Distributed  Active  Archive  Center  (LP  DAAC)  at 
the  USGS/Earth  Resources  Observation  and  Science  (EROS)  Center,  Sioux  Falls, 
South  Dakota.  The  data  product  for  the  image  was  provided  by  NASA  [12]. 

2.3.3  Cloud-Aerosol  LIDAR  and  Infrared  Pathfinder  Satellite  Observations 
(CALIPSO) 

The  CALIPSO  satellite  is  also  a  member  of  the  NASA  A-Train,  following  about 
one  minute  and  15  seconds  behind  Aqua.  CALIPSO ’s  LIDAR  emits  and  measures 
backscatter  at  532  nm  and  1064  nm.  This  instrument  has  the  ability  to  not  only  detect  the 
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top  and  bottom  heights  of  multiple  cloud  layers  but  also  to  distinguish  between  aerosol 
particles  and  cloud  droplets  and  retrieve  information  about  cloud  properties.  The  ability  to 
differentiate  between  effects  of  aerosols  and  clouds,  and  the  determination  of  cloud  type, 
allows  accurate  estimation  of  wavelength  dependent  path  extinction.  Figure  7  shows  an 
example  of  how  CALIPSO’s  vertical  resolution  and  cloud  property  measurements  can  be 
combined  with  data  from  other  instruments.  In  this  case  the  2-D  CALIPSO  data  are 
overlaid  on  a  corresponding  Moderate  Resolution  Imaging  Spectroradiometer  (MODIS) 
image,  both  collected  as  the  A-Train  passed  over  Hurricane  Bill  in  2009. 


Figure  7.  NASA  image  of  overlapping  CALIPSO  and  MODIS  data— ascending 
satellite  pass  over  Hurricane  Bill  in  2009  (Courtesy  NASA)  [13]. 
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2.4  Numerical  Weather  Prediction  (NWP) 


2.4.1  Global  Forecast  System  (GFS) 

GFS  is  a  global,  hydrostatic  weather  forecast  model  produced  by  the  National 
Centers  for  Environmental  Prediction  (NCEP).  The  current  GFS  Global  Spectral  Model  is 
operated  at  its  standard  0.5  degree  horizontal  resolution  grid,  but  an  improved  0.25  degree 
grid  became  operational  in  late  2014.  Both  models  are  run  4  times  per  day,  each  with  3 
hour  forecast  timesteps  out  to  240  hours,  followed  by  a  lower  temporal  resolution  extension 
period  at  12  hour  timesteps  out  to  384  hours.  The  model  runs  on  64  vertical,  unequally- 
spaced  sigma  levels.  This  means  the  pressure  levels  in  the  model  are  normalized  by 
dividing  level  pressure  by  the  surface  pressure,  cr  =  p /  ps ,  compressing  the  levels  so  that 

none  of  the  model  pressure  levels  exist  below  the  Earth’s  surface.  The  64  model  levels 
extend  up  to  a  highest  pressure  level  of  0.03  hPa,  but  the  output  product  only  includes  26 
levels  with  a  top  pressure  level  of  10  hPa.  For  this  reason,  the  top  of  the  GFS  coverage 
volume  in  the  available  data  sets  is  typically  around  30-32  km.  The  gridded  model  output 
is  publically  available  for  download  through  the  NOAA  Operational  Model  Archive  and 
Distribution  System  (NOMADS)  server  [14]  [15]. 

2.4.2  European  Centre  for  Medium-Range  Weather  Forecasts  (ECMWF) 

An  NWP  model  that  provides  a  source  for  temperatures,  heights  of  pressure  levels, 
and  winds  above  the  top  of  the  available  GFS  data  is  the  European  Centre  for  Medium- 
Range  Weather  Forecasts  (ECMWF).  Of  particular  interest  in  this  research  is  the  ECMWF 
Re-Analysis  (ERA)  project  ERA-20C.  The  ERA-20C  data  set  represents  output  from  a 
global,  coupled  atmosphere,  land  surface,  and  ocean  wave  model,  which  was  re-run  with 
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assimilated  surface  observations.  The  data  set  has  a  3-hour  temporal  resolution  and  91 
vertical  model  levels  up  to  a  top  level  pressure  of  0.01  hPa,  above  75  km,  covering  the 
entire  troposphere,  stratosphere,  and  much  of  the  mesosphere.  The  ECMWF  re-analysis 
data  are  available  to  the  public,  covering  dates  from  the  year  1900  through  the  end  of  2010, 
overlapping  the  first  eight  years  of  AIRS  data  [16]. 

This  ECMWF  model  is  non-hydrostatic,  meaning  that  it  is  able  to  calculate  vertical 
pressure  gradients  and  accelerations.  Hydrostatic  models  assume  hydrostatic  equilibrium, 
forcing  the  downward  weight  of  the  atmosphere  to  balance  the  upward  pressure  gradient 

force.  At  all  points  in  hydrostatic  models,  =  ~pg  ,  where  p  is  pressure,  z  is  height,  p 

is  density  of  the  air,  and  g  is  acceleration  due  to  gravity  at  the  height  z.  Large  scale  vertical 
motion  can  be  deduced  in  hydrostatic  models  from  the  convergence  or  divergence  of  air 
flows,  but  pressure  perturbations  not  in  balance  with  gravity  are  not  allowed.  Non¬ 
hydrostatic  models  include  additional  equations,  not  present  in  hydrostatic  models.  These 
extra  equations  allow  smaller  scale  vertical  accelerations  to  be  directly  calculated,  but  their 
inclusion  comes  at  the  cost  of  making  the  model  more  computationally  expensive  to  run. 
The  vertical  acceleration  equations  are  necessary  for  the  model  to  be  able  to  accurately 
account  for  buoyancy  forces,  convection,  and  any  local  deviations  from  hydrostatic 
equilibrium  [15]  [17] 

The  global  coverage  and  excellent  horizontal  and  vertical  resolution  of  the  ECMWF 
data  make  them  valuable  for  AIRS  data  comparison  and  Cn2  calculation.  Figure  8,  Figure 
9,  and  Figure  10  show  sample  ECMWF  output  for  1800  UTC  on  10  Mar  2010. 


18 


ECMWF  Temperature  (K)  803  hPa  Level 


Figure  8.  ECMWF  Temperature,  803  hPa  level,  10  Mar  2010, 1800  UTC. 


ECMWF  Cloud  Fraction  803  hPa  Level 


Longitude 


Figure  9.  ECMWF  Cloud  Fraction,  803  hPa  level,  10  Mar  2010, 1800  UTC. 
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Vertical  Wind  Component  (m/s)  803  hPa  Level 


Figure  10.  ECMWF  vertical  wind  component,  803  hPa  level,  10  Mar  2010, 1800 
UTC. 

2.4.3  Comparison  of  NWP  models  investigated 

The  primary  NWP  sources  utilized  in  this  research  are  the  GFS  and  ECMWF 
models.  Additional  models  were  studied  but  were  not  included  in  the  Cn2  derivation 
methods  developed.  The  ability  to  run  the  Weather  Research  and  Forecasting  (WRF) 
model  internally  (with  a  customizable  domain  and  grid)  and  the  1  hour  timestep  of  the  High 
Resolution  Rapid  Refresh  (HRRR)  model  are  both  advantages,  but  the  fact  that  neither  of 
these  models  provide  global  coverage  led  to  them  not  being  incorporated  into  the 
techniques  developed.  Table  1  details  the  coverage  and  resolution  of  five  of  the  NWP 
models  investigated  for  use  in  this  research. 
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Table  1.  Resolution  comparison  for  NWP  models  investigated 


Model 

Domain 

Horizontal 

Grid 

Spacing 

Vertical 

Levels 

Top  Level 

Pressure 

Temporal 

Resolution 

Non- 

Hydrostatic 

Global  Forecast  System  (GFS)  0.50  deg 

Global 

0.50  deg 

26 

10hPa 

3  hr 

No 

Global  Forecast  System  (GFS)  0.25  deg 

Global 

0.25  deg 

26 

lOhPa 

3  hr 

No 

Weather  Research  and  Forecasting  (WRF) 

Customizable  to  meet  requirements 

Yes 

High  Resolution  Rapid  Refresh  (HRRR) 

CONUS 

3  km 

50 

20  hPa 

1  hr 

Yes 

ECMWF  ERA-20C  Reanalysis 

Global 

0.125  deg 

91 

0.01  hPa 

3  hr 

Yes 

2.5  Atmospheric  Turbulence  Theory 

Optical  turbulence  in  the  atmosphere  poses  a  problem  for  most  electro-optical 
systems.  It  leads  to  image  distortion  for  imaging  sensors  and  beam  spread  for  laser 
systems.  Random  fluctuations  in  the  refractive  index  of  air  cause  optical  turbulence, 
leading  to  variations  in  the  speed  at  which  a  wavefront  propagates.  The  quantitative 
measure  of  the  strength  of  the  optical  turbulence  that  is  most  often  used  is  the  refractive 
index  structure  function  Cn2  [18]. 

A  structure  function  describes  the  spatial  statistics  of  the  parameter  of  interest  and 
can  be  defined  as 

5>0)  =  ([<zK*)-(zK-*:')]  )  (6) 


where  the  separation  between  the  two  positions  is  designated  r, 

r  =  \x  —  x'\ 


(7) 


The  brackets  (  }  in  Equation  (6)  denote  an  ensemble  average  and  isotropy,  and 
homogeneity  of  the  turbulence  variance  are  implied  [19]. 
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Within  the  inertial  subrange  (the  range  of  scales  where  turbulence  is  neither 
generated  nor  dissipated,  bounded  by  the  inner  scale  lo  and  the  outer  scale  Lo)  Kolmogorov 
turbulence  follows  the  “two-thirds  power  law”  [20].  For  two  positions  separated  by  a 
distance  r,  the  relationship  between  the  index  of  refraction  n  at  each  position  and  the  index 
of  refraction  structure  function  C„ 2  can  be  expressed  as 

|[fi(x)-R(x  +  r)]2|  =  CV/3.  (8) 

In  this  inertial  subrange,  C2  is  also  proportional  to  the  level  of  the  power  spectrum  of 
refractive  index  fluctuations 

<D„(fci)  =  0.25  C„V/3  (9) 

where  ®n(^)  is  the  average  one-dimensional  power  spectrum  of  the  index  of  refraction  n, 
and  kt  is  the  wave-number  component  in  the  direction  of  propagation  [21]  [22],  The 

structure  function  for  refractive  index,  Cn2,  provides  a  measure  of  the  strength  of  the  optical 
turbulence.  Typical  Cn2  values  for  strong  optical  turbulence  exceeds  10"13  m'2/3  and  weak 
turbulence  has  C2  values  smaller  than  10'16  m"2/3  [23]  [24]  [25]. 

Atmospheric  flow  is  complex  in  part  because  turbulent  eddies  can  exist 
simultaneously  on  all  spatial  scales  from  the  planetary  scale  (on  order  of  20,000  km)  down 
to  the  fine  scale  (on  order  of  2  mm)  and  viscous  dissipation  subrange  (molecular  scale). 
The  resulting  atmospheric  flow  is  the  superposition  of  all  of  these  turbulent  eddies  [26]. 
There  is  also  complex  interaction  between  horizontal  motion  and  vertical  motion. 
Buoyancy  effects  caused  by  heat  transfer  or  phase  change  can  have  an  effect  on  the  type 
and  strength  of  the  turbulence  [21].  The  energy  source  for  all  of  this  motion  is  the  non- 
uniform  global  distribution  of  differential  heating,  and  the  turbulent  energy  is  transferred 
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from  each  scale  to  the  next  smaller  scale,  from  planetary  scale  motion  like  the  jet  stream 
down  to  the  smallest  eddies  where  it  is  lost  to  friction  [27].  The  three  spectral  regions 
commonly  discussed  in  atmospheric  turbulence  research  are  [21]  [24]  [28]: 

1)  Energy  production  range — turbulent  energy  is  produced  by  buoyant  forces  and 
mechanical  shear  (on  order  of  2  km). 

2)  Inertial  subrange — no  energy  is  produced  or  dissipated;  it  is  transferred  from 
large  turbulent  eddies  to  smaller  eddies  (scales  between  the  outer  scale  Lo  and 
the  inner  scale  lo). 

3)  Dissipation  range — turbulent  kinetic  energy  is  converted  to  heat  through  friction 
(on  order  of  1  cm  and  smaller). 

The  work  of  Kolmogorov  provides  a  unifying  view  of  turbulence  dynamics  across 
these  spectral  regions.  Kolmogorov’s  theory  describes  small  scale  turbulence  structure  as 
statistically  homogeneous  and  isotropic  [23].  The  distribution  of  turbulent  energy  within 
the  inertial  subrange,  according  to  Kolmogorov’s  hypothesis,  can  be  determined  by  the  5/3 
law,  which  states  subrange  spectral  densities  are  proportional  to  k~5/3 ,  where  k  is  the 
magnitude  of  the  three-dimensional  wave  number.  These  spectral  densities  allow  the 
determination  of  the  energy  transfer  rate  or  energy  dissipation  rate  [24]  [29]. 

While  Cn2  values  vary  significantly  with  changes  in  atmospheric  conditions,  there  is 
a  general  dependence  on  height  that  allows  for  the  creation  of  standard  vertical  profile 
models.  One  of  these  standard  models  is  Hufnagel- Valley,  with  the  common 
implementation  H-V  5/7.  In  this  implementation,  when  integrated  for  a  vertical  path  at  a 
wavelength  of  /l=0.5  pm,  the  resulting  turbulent  coherence  length  po  is  5  cm,  and  the 
isoplanatic  patch  size,  6o,  is  7  prad.  These  values  are  typical  for  nighttime  observations  in 
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the  field  of  astronomy.  Turbulent  coherence  length,  the  length  over  which  phase  distortion 
has  the  mean  square  value  of  1  rad2,  is  a  parameter  used  to  describe  the  amount  of  blur 
effects  on  receiver  systems.  The  parameter  isoplanatic  patch  size  also  comes  from  the  field 
of  astronomy,  where  it  defines  the  maximum  angular  distance  at  which  two  moving  points 
appear  to  shift  position  in  unison.  At  this  angular  difference  between  paths,  the  mean 
square  wavefront  error  is  1  rad2.  [30]  [31].  H-V  has  the  form 

C„2(z)  =  8.148xl0“26£/2zxVz  +2.700xl0“16e“z/1'5  +C„2(OK10z  (10) 

where  the  surface  value  Cn2( 0)  =  1.7xl0"14  m"2/3  and  an  upper  level  wind  speed  U  =  21  m/s 
yield  the  5/7  implementation  [28]  [30]  [32]  [33],  Figure  11  depicts  H-V  5/7  plotted  against 
another  common  standard  model,  Critical  Laser  Enhancing  Atmospheric  Research 
(CLEAR)  1  Night. 


Figure  11.  Comparison  of  two  standard  vertical  Cn2  profile  models,  H-V  5/7  and 
CLEAR  1  Night. 
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The  refractive  index  of  air  is  a  function  of  pressure,  temperature,  and  water  vapor 
concentration.  At  radio  wavelengths  the  most  important  of  these  factors  is  the  water  vapor. 
In  analysis  of  refractive  index  fluctuations  affecting  radio  transmission,  it  is  common  to 
consider  only  the  water  vapor  concentration  gradients  or,  in  the  approach  described  by 
Tatarskii,  to  consider  the  water  vapor  and  temperature  gradients.  At  optical  wavelengths 
the  water  vapor  has  little  influence  on  refractive  index,  and  the  standard  approach  taken  is 
to  base  optical  refractive  index  fluctuations  solely  on  temperature  fluctuations.  In  both  of 
these  regions  of  the  electromagnetic  spectrum,  the  traditional  approach  is  to  treat  pressure 
fluctuations  as  negligible  [34]  [35].  A  pressure  gradient  term  is  derived,  analyzed,  and 
incorporated  into  a  modified  Cn2  calculation  process  detailed  in  Section  4.3.5. 
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3  Methodology 


3.1  AIRS  Temperature  Measurements  and  Retrieving  AIRS  Data 

AIRS  data  files  are  made  available  to  the  public  by  NASA,  searchable  by  time  and 
geographic  coverage  through  the  Mirador  Earth  science  data  search  tool  developed  and 
maintained  by  GES-DISC  [36].  For  this  research,  an  automated  process  built  using 
MATLAB®  is  used  to  search  the  NASA  database  and  retrieve  the  required  AIRS  data  files 
for  processing.  In  order  to  have  radiosonde  wind  profiles  available  to  compare  to  the 
profiles  derived  from  satellite  soundings,  many  of  the  locations  for  this  research  were 
selected  in  the  vicinity  of  rawinsonde  sites  that  the  Aqua  satellite  passes  over  near  the 
scheduled  daily  weather  balloon  release  times  (0000  UTC  and  1200  UTC)  [37].  These 
locations  lie  in  north-south  bands  near  20°  E  and  160°  W  longitude,  and  due  to  Aqua’s  sun- 
synchronous  orbit,  the  location  of  these  bands  remains  constant  throughout  the  year. 
Rawinsonde  observations  for  this  research  were  obtained  through  a  website  maintained  by 
the  University  of  Wyoming’s  Department  of  Atmospheric  Science  [38]. 

The  continuous  stream  of  AIRS  data  is  broken  into  six-minute  granules.  Each 
AIRS  data  point  within  the  granule  corresponds  to  one  of  the  1350  sounding  locations 
covered  during  those  six  minutes  (45  points  along  the  flight  path,  30  points  wide).  The 
field-of-view  center  coordinates  for  each  of  the  sounding  locations  in  one  data  granule  are 
shown  in  Figure  12.  This  file  covers  a  segment  of  a  descending  AIRS  pass  that  tracked 
along  the  eastern  border  of  Alaska. 
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AIRS  Field-of-View  Boresight  Coordinates 
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Figure  12.  Depiction  of  sounding  boresight  coordinates  for  a  6-minute  data  granule. 
In  this  example  granule  number  125  on  15  Oct  2014  represents  this  segment  of  this 
descending  AIRS  pass  over  Alaska  between  1229  and  1235  UTC. 

The  level  2  AIRS  data  set  contains  atmospheric  temperature  values  for  100  standard 
vertical  pressure  levels  at  each  sounding  location.  These  pressure  levels  are  the  same  for 
each  sounding  column,  making  it  straightforward  to  compare  one  sounding  to  another.  The 
first  step  in  selecting  the  appropriate  AIRS  files  for  download  is  the  determination  of  a  date 
and  time,  and  a  set  of  coordinates  of  interest.  The  Mirador  server  is  searched,  and  if  there 
is  an  AIRS  granule  that  includes  the  entered  coordinates  during  the  specified  time  search 
window,  it  is  downloaded  (two  granules  are  downloaded  if  the  entered  coordinates  fall 
within  300  km  of  the  leading  or  trailing  edge  of  the  granule).  For  one  of  these  data  granules, 
a  horizontal  layer  of  temperatures  at  the  AIRS  standard  pressure  level  of  753.6  hPa  is 
shown  in  Figure  13. 
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Figure  13.  Horizontal  slice  through  the  volume  of  AIRS-retrieved  temperatures 
contained  in  a  6-minute  data  granule  at  753.6  hPa,  1229-1235  UTC,  15  Oct  2014. 

The  retrieved  data  are  filtered  to  select  the  vertical  temperature  profiles  for  all 
sounding  locations  within  300  km  of  the  coordinates  of  interest.  This  filtered  data  set 
includes  approximately  160-200  sounding  locations  within  each  pressure  layer,  making  up 
a  circular  region  600  km  across.  For  each  sounding  location,  the  data  set  includes  surface 
pressure,  surface  geopotential  height,  and  a  temperature  measurement  for  up  to  100 
pressure  levels  between  0.0161  hPa  and  1100  hPa  (all  100  temperatures  may  not  be 
available  due  to  terrain  elevation  and  surface  pressure).  During  the  level  2  processing  to 
calculate  temperature  values  from  the  satellite  measured  IR  radiances  and  microwave 
brightness  temperatures,  an  error  value  is  assigned  to  each  temperature. 

The  AIRS  field-of-view  coordinates  provide  the  correct  latitude  and  longitude 
coordinates  for  only  the  lowest  altitude  field-of-view,  where  the  AIRS  instrument’s  line  of 
sight  intersects  the  surface  of  the  Earth.  As  the  instrument  sweeps  across  each  scanline, 
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the  angle  varies  between  49  degrees  left  of  nadir  to  49  degrees  right  of  nadir.  The 
coordinates  of  each  temperature  measurement  above  the  surface  layer  are  not  calculated, 
nor  are  they  provided  with  the  AIRS  level  2  products,  but  the  satellite’s  latitude,  longitude, 
and  altitude  are  reported  throughout  each  orbit.  The  atmospheric  height  of  each 
temperature  measurement  is  known,  so  extending  a  line  from  each  surface  location  to  the 
satellite’s  position  allows  assignment  of  correct  latitude  and  longitude  coordinate  to  each 
measurement  along  that  line.  Figure  14  depicts  the  result  of  this  process,  showing  the 
corrected  location  of  each  of  the  135,000  temperature  retrievals  during  a  6-minute  data 
granule.  Each  red  line  in  this  figure  depicts  the  column  of  100  soundings,  which 
correspond  to  a  single  AIRS  field-of-view  at  the  surface. 


Figure  14.  Corrected  latitude  and  longitude  coordinates  for  entire  3-D  granule  of 
AIRS  temperature  data  during  an  ascending  AIRS  pass  over  Florida. 


29 


3.2  AIRS-Derived  Wind  Profiles — Vertical  Wind  Profile  Calculation  Using  the 
Thermal  Wind  Equations 

For  each  of  the  100  AIRS  pressure  levels,  smooth  surfaces  are  fit  through  the 
pressure  heights  and  temperature  distributions.  These  surfaces  are  used  to  find  the  local 
slope  of  the  isobaric  surface  and  the  rate  of  temperature  change  along  the  pressure  level  in 
both  the  east-west  and  the  north-south  directions.  An  example  AIRS  data  granule  is  plotted 
in  Figure  15. 
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Figure  15.  Sample  of  sounding  locations  surrounding  coordinates  of  interest — along 
satellite  flight  path  (blue  diamonds)  and  perpendicular  to  satellite  flight  path  (red 
squares).  This  granule  of  AIRS  data  was  collected  on  24  Oct  2014,  as  the  Aqua 
satellite  overflew  Dayton,  OH  at  1839  UTC.  This  horizontal  slice  through  the  data 
shows  retrieved  temperature  at  the  151  hPa  pressure  level  (heights  between  13.1  and 
14.2  km). 
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Figure  16  and  Figure  17  illustrate,  in  two  dimensions,  how  the  smoothing  process 
leads  to  the  required  gradients  [1].  Figure  16  shows  two  cross-sections  through  a  sounding 
data  set.  The  north-south  trail  of  blue  data  points  show  a  slice  of  sounding  locations  parallel 
to  the  flight  path  of  the  satellite  and  the  red  locations  identify  an  east-west  scanline 
perpendicular  to  satellite  track,  corresponding  to  the  blue  and  red  tracks  in  Figure  15.  The 
two  tracks  intersect  at  the  sounding  location  nearest  to  the  input  target  coordinates.  In 
Figure  16,  the  temperature  measurements  within  the  151  hPa  layer  for  each  sounding 
location  along  these  two  lines  are  plotted.  The  small  error  and  variation  in  the  temperatures 
can  lead  to  many  different  values  for  temperature  gradient,  depending  on  how  many  data 


points  are  included,  and  how  a  surface  is  fit  to  the  data. 
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Figure  16.  Satellite  measured  temperatures  at  the  151  hPa  pressure  level  along  the 
satellite  flight  path  (blue  diamonds)  and  perpendicular  to  the  satellite  flight  path 
(red  squares)  for  the  data  granule  shown  in  Figure  15.  The  vertical  dotted  lines 
indicate  the  coordinates  of  interest. 
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To  smooth  the  data,  and  calculate  the  most  accurate  temperature  and  height 
gradients,  a  surface  needs  to  be  fit  to  the  data  to  reduce  the  effect  of  error  in  the  temperature 
and  height  data.  A  three-dimensional  surface  is  fit  to  each  level’s  data  by  local  regression 
smoothing  (specifically  the  Robust  Locally  Weighted  Scatter  Plot  Smoothing  function 
within  the  MATLAB®  Curve  Fitting  Toolbox  [39]).  Representative  cross-sections  of  this 
smooth  surface  are  shown  in  Figure  17.  The  temperature  gradient  is  determined  by  finding 
the  slope  of  this  surface  at  the  target  coordinates,  and  by  a  similar  curve  fitting  process,  the 
slopes  of  the  isobaric  surfaces  are  determined.  This  method  of  smoothing  accounts  for 
more  data  points  than  just  the  few  immediately  surrounding  the  target  coordinates,  and  it 


minimizes  the  effect  of  error  in  the  individual  satellite  measured  temperatures. 
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Figure  17.  Cross-sections  of  the  3-D  smoothed  surface  fit  to  all  AIRS  temperature 
measurements  within  300  km  of  the  input  target  coordinates  for  the  151  hPa 
pressure  level  data  shown  in  Figure  15  and  Figure  16. 
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The  vertical  temperature  gradients  are  found  by  dividing  the  difference  in 
temperature  at  the  center  of  two  pressure  surfaces  by  the  vertical  separation  between  the 
surfaces,  and  repeating  this  process  for  each  pair  of  pressure  surfaces.  Once  all  of  the 
required  temperature  and  pressure  height  gradients  are  calculated,  the  wind  profile  can  be 
assembled  from  the  bottom  up  or  the  top  down,  one  level  at  a  time. 

The  thermal  wind  relationship  between  components  of  the  geostrophic  wind,  ug  and 
vg,  at  two  consecutive  levels  in  the  atmosphere,  1  and  2,  is  described  in  Equations  (11)  and 
(12).  This  development  makes  two  important  assumptions — that  the  atmosphere  is  in  both 
hydrostatic  and  geostrophic  balance.  When  expressed  as  components,  u  is  the  zonal 
velocity,  the  component  of  the  horizontal  wind  from  the  west,  and  v  is  the  meridional 
velocity,  the  component  of  the  horizontal  wind  from  the  south.  The  change  in  magnitude 
of  each  wind  component,  from  one  level  to  the  next,  can  be  found  by 

(“T -(“«),=-  (11) 

and 


(v  )  -tv) 

[  sj2  \  sJi  f  dx 


(12) 


where  go  is  the  acceleration  due  to  gravity  at  sea  level,  and  /  is  the  Coriolis  parameter, 
which  depends  on  latitude.  The  thickness  of  each  layer  z2  ~zt  is  linearly  proportional  to 


the  mean  temperature  within  the  layer,  as  described  by  the  hypsometric  equation  [26]: 
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where  Rd  is  the  dry  gas  constant,  287  J/K  kg,  and  Tv  is  the  mean  virtual  temperature  in  the 
layer.  The  variables  pi  and  pi  are  the  atmospheric  pressures  at  the  levels  at  the  bottom  and 
top  of  the  layer.  Combining  Equations  (11),  (12),  and  (13)  provides  a  method  for 
determining  the  vertical  change  in  u  and  v  components  of  the  geostrophic  wind  through 
each  layer.  The  change  in  wind  between  the  top  and  bottom  of  a  pressure  layer  is  a  function 
of  the  vertical  change  in  pressure  through  the  layer  and  the  horizontal  temperature  gradient, 
as  described  in  Equations  (14)  and  (15): 
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The  form  of  the  thermal  wind  relationship  shown  in  Equations  (14)  and  (15) 
neglects  the  dependence  on  vertical  temperature  gradient,  which  also  contributes  to  the 
wind  gradient  [26].  A  more  complete  set  of  equations  for  the  wind  gradients  in  u  and  v 
components  is  [40]: 
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Starting  from  one  known  wind  level,  in  this  case  the  GFS  model  surface  wind,  the 
full  wind  profile  is  calculated  by  applying  Equations  (16)  and  (17)  to  each  layer  in 
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sequence,  adding  the  change  in  wind  component  through  the  layer  to  the  value  at  the 
previous  level.  This  technique  performs  better  when  calculating  total  wind  speed  than  wind 
direction.  Surface  friction  effects  can  cause  rotation  of  the  wind  vector  that  is  not  reflected 
in  the  temperature  data.  Analysis  of  how  consistent  and  accurate  the  AIRS -derived  wind 
profiles  are  is  presented  in  Section  4.2.1.  Figure  18  and  Figure  19  show  example  AIRS- 
derived  wind  profiles,  compared  to  the  corresponding  rawinsonde  observations  (RAOB) 
profiles. 


Figure  18.  Comparison  of  AIRS-derived  and  RAOB  u  and  v  wind  components  for 
11  Apr  2014  1250  UTC,  Anchorage,  AK. 
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Figure  19.  Comparison  of  AIRS-derived  and  RAOB  wind  speed  and  direction  for 
11  Apr  2014  1250  UTC,  Anchorage,  AK  (same  wind  profile  as  Figure  18). 

3.3  Calculation  of  the  Index  of  Refraction  Structure  Function  Cn2 


3.3.1  Relationship  between  Ct2  and  C,,2 
The  AIRS  instrument  provides  the  3-D  temperature  data  needed  for  calculation  of 
the  temperature  structure  functions  Ct2  and  C,,2  for  optical  wavelengths.  Using  an 
expression  derived  by  Tatarskii  [35]  [41]  [42],  Ct  can  be  determined  by 

f  dT 


CT  =  a  aL 07 


dz 


+  rd 


(18) 


where  a2  is  an  empirically  determined  constant,  most  often  assigned  the  value  of  2.8  [43] 
[44]  [45]  [46],  The  outer  scale  length  L0  is  the  largest  scale  of  the  atmospheric  turbulence, 
or  the  upper  bound  of  the  inertial  subrange.  Values  for  L0  are  reported  throughout  the  range 
of  “many  meters  to  several  tens  of  meters,”  but  experimental  measurements  are  scarce, 
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especially  outside  the  planetary  boundary  layer  [47].  Unless  noted  otherwise,  L0  is 
estimated  to  be  100  m  in  this  research.  The  dry  adiabatic  lapse  rate  yd  is 


yd  =^-  =  9.8x10  ~3K/m 

cp 


(19) 


where  g  is  acceleration  due  to  gravity  at  sea  level,  9.81  m/s,  and  cp  is  the  heat  capacity  of 

air  at  constant  pressure,  1005  J/kg  K  for  dry  air  [5],  The  term  +  yd  j  is  the  deviation 

between  the  temperature  lapse  rate  and  the  dry  adiabatic  lapse  rate,  making  it  another  way 
to  express  the  vertical  rate  of  change  of  potential  temperature  6.  Potential  temperature  is 
the  temperature  a  parcel  of  unsaturated  air  would  acquire  if  it  were  brought  adiabatically 
to  a  standard  pressure,  usually  1000  hPa.  Potential  temperature  can  be  calculated  by 


d  =  T 
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where  p0  is  the  standard  pressure,  1000  hPa,  and  the  ratio  of  the  dry  gas  constant  to  the  heat 
capacity  of  dry  air,  Rd/cp  ~  0.286  [5], 

The  coefficient  a  in  Equation  (18)  can  be  replaced  by  the  ratio  of  eddy  coefficients 
Kh  and  Km  by  a  «  KH  / KM  .  KH  is  the  eddy  diffusivity  for  heat,  and  KM  is  the  eddy 

diffusivity  for  momentum.  The  ratio  of  these  eddy  transfer  coefficients,  Kh/Km  ,  is  also 

referred  to  as  the  inverse  of  the  turbulent  Prandtl  number.  By  replacing  a  with  this  ratio, 
Equation  (18)  becomes 
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Cn2  can  be  related  to  Cr  by  [28]  [35]  [44] 

( p„  \2  (  p\2 

C2  =  —  Cl=  79xl0-6  — y  C2  .  (22) 

"  ydT )  T  y  T2) 

The  numerical  coefficient  79xl0“6  in  Equation  (22)  is  appropriate  for  pressure  measured 
in  hPa  and  temperatures  in  K.  This  coefficient  is  wavelength  dependent,  and  the  value 
79x10  6  is  commonly  selected  because  it  corresponds  to  the  middle  of  the  visible  optical 
range.  The  value  would  be  82.9  xio-6  at  300  nm  and  77.4  xio-6  at  700  nm  [35]. 

Combining  Equations  (21)  and  (22)  yields  the  initial  method  for  calculating  Cn2 
employed  in  this  research,  which  will  be  referred  to  in  this  document  as  the  traditional 
Tatarskii  approach: 

K  f  p(2  4//'  pt  \2 

C„2=2.8-^  79x10-6^  W3  ^~  +  rd  ■  (23) 

AJlf  v  1  J  \OZ  j 

3.3.2  Gradient  Richardson  Number 

The  gradient  Richardson  number,  Ri,  indicates  the  dynamic  stability  of  the 
atmosphere  and  gives  insight  into  turbulence  production.  It  is  the  ratio  of  buoyancy  forces 
(the  numerator  is  the  square  of  the  Brunt- Vaisala  frequency)  to  inertial  accelerations  (the 
denominator  is  the  square  of  the  vertical  shear  in  the  horizontal  wind).  It  can  be  calculated 
from  the  vertical  gradients  in  potential  temperature  and  wind  speed,  previously  derived 
from  AIRS  temperature  data  using  [26]  [34]  [48]  [49] 
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In  Equation  (24)  Ov  is  the  virtual  potential  temperature,  z  is  the  geometric  height,  and  V  is 
the  magnitude  of  the  horizontal  wind.  The  Richardson  number  is  a  measure  of  the  stability 
of  the  flow,  with  larger  values  corresponding  to  stable,  laminar  flow.  When  the  Richardson 
number  exceeds  a  critical  value,  usually  reported  to  be  near  0.2  or  0.3,  turbulent  flow 
becomes  laminar  [21]  [47]  [50],  Table  2,  reproduced  from  [34],  identifies  typical  gradient 
Richardson  numbers  for  various  conditions  of  static  stability  and  flow  type. 


Table  2.  Values  of  the  gradient  Richardson  number  for  different  static  stability  and  flow 
_ conditions,  reproduced  from  Avila  and  Vernin  [34]. _ 
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turbulent 
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<0 
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The  eddy  diffusivity  ratio,  Kh/Km  ,  can  be  approximated  from  the  gradient 
Richardson  number  by  [51]: 
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3.3.3  Application  of  the  Traditional  Tatarskii  Approach  to  AIRS  Data 

This  Cn2  calculation  process  has  been  previously  applied  to  Numerical  Weather 
Prediction  (NWP)  model  data.  Alliss  and  Felton  used  the  WRF  model  to  generate  Cn2 
climatology  over  the  Hawaiian  Islands  in  2009  [44]  [45],  but  the  techniques  developed  to 
apply  the  process  to  satellite  sounding  data,  where  the  required  temperature,  pressure,  and 
wind  information  is  all  passively  collected  from  space,  has  not  been  found  in  the  literature 
on  the  subject.  An  example  of  a  Cn  profile  calculated  by  applying  this  process  to  AIRS 
temperatures  and  derived  winds  is  shown  in  Figure  20.  Additional  results  and  analysis  of 
the  technique  are  found  in  Section  4.3. 


Vertical  C2  Profile 


Turbulence  Structure  C2  (m'2/3) 

Figure  20.  AIRS-derived  Cn2  profile  - 11  Apr  2014  1250  UTC,  Anchorage,  AK. 
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3.4  Cloud  Location  and  Characterization 


3.4.1  Determination  of  Cloud-Top  Height  using  AIRS  Data 

An  additional  focus  of  this  research  is  determination  of  the  heights  of  cloud  layers  and 
tops.  Clouds  cover  approximately  40%  of  the  Earth’s  surface  at  any  given  time  [52],  At 
wavelengths  of  interest  for  directed  energy  applications,  direct  transmittance  through  all 
but  the  thinnest  of  cloud  layers  is  extremely  low  [2],  For  this  reason,  it  is  essential  to  be 
able  to  identify  whether  a  specified  slant-range  path  through  the  atmosphere  is  cloud-free 
or  not.  A  critical  step  in  this  direction  is  to  identify  the  height  of  the  top  of  the  highest 
cloud  layer  at  each  sounding  location,  which  is  related  to  the  greatest  atmospheric  pressure 
for  which  AIRS  is  able  to  retrieve  a  temperature.  Satellite-based  sounders  provide  the  data 
necessary  to  both  locate  the  extent  of  the  cloud  cover  horizontally  and  determine  the  height 
of  the  highest  layer.  Thermal  radiation  emitted  from  heights  below  clouds  is  largely 
absorbed  before  reaching  space,  and  the  channels  whose  weighting  functions  are  associated 
with  those  heights  measure  unrealistically  cold  temperatures.  Figure  21  depicts  an  image 
of  the  cloud  cover  over  the  eastern  United  States  from  one  of  the  Geostationary  Operational 
Environmental  Satellites  (GOES),  and  Figure  22  shows  plots  of  the  AIRS  cloud  fraction 
and  cloud  top  pressure  for  one  6-minute  granule  of  data  from  an  ascending  pass  about  10 
minutes  after  the  GOES  data  were  collected.  The  left  plot  shows  the  AIRS-derived  cloud 
fraction,  the  fraction  of  the  nine  included  AIRS  fields-of-view  that  are  not  cloud-free  all 
the  way  to  the  surface.  The  right  plot  shows  the  atmospheric  pressure  for  the  AIRS 
standard  layer  containing  the  highest  layer  determined  to  be  cloudy,  allowing  an  estimate 
of  the  altitude  of  the  top  of  the  highest  cloud  layer. 
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Figure  21.  GOES-East  image  for  1815  UTC  on  27  Jan  2014,  retrieved  from  the 
NOAA  Geostationary  Satellite  Server  [53]. 
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Figure  22.  AIRS-derived  cloud  fraction  (left)  and  cloud-top  pressure  in  hPa  (right). 
AIRS  data  collected  between  1824  and  1830  UTC,  27  Jan  2014  [41]. 
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A  3-D  depiction  of  the  cloud-top  surface  that  results  from  assigning  heights  to  the 
AIRS-derived  cloud-top  pressures  is  shown  in  Figure  23.  Holes  or  missing  data  in  the 
surface  indicate  regions  that  are  cloud-free  all  the  way  to  the  surface  [54], 
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Figure  23.  AIRS-derived  cloud-top  heights  - 1826-1830  UTC,  27  Jan  2014. 

3.4.2  Determination  of  Cloud-Top  Height  using  MODIS  Data 

By  incorporating  data  from  a  second  A-Train  instrument,  cloud  layer  identification 
and  top  height  calculation  can  achieve  a  level  of  accuracy  and  resolution  that  will  make  the 
near  real-time  cloud  characterization  useful  for  directed  energy  research  and  operations. 
MODIS  data  is  available  from  the  instruments  onboard  the  Aqua  and  Terra  satellites  [55], 
The  1  km  horizontal  resolution  of  the  MODIS  instrument  can  improve  the  determination 
of  the  horizontal  extent  of  the  cloud  and  better  identify  locations  that  are  cloud-free  all  the 
way  to  the  surface.  Figure  24  shows  the  MODIS  cloud  products  for  the  same  geographic 
area  and  time  period  covered  by  the  AIRS  data  in  Figure  22  and  Figure  23. 
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Figure  24.  MODIS-derived  cloud-top  pressure  in  millibars  (left)  and  cloud-top 
height  in  meters  (right).  Collected  between  1825-1830  UTC,  27  Jan  2014. 

Figure  25  shows  a  MODIS-derived  3-D  depiction  of  cloud-top  heights,  similar  to 
the  AIRS-derived  image  shown  in  Figure  23.  The  improvement  in  horizontal  resolution  is 
clear,  with  the  1-km  resolution  better  defining  cloud-top  shape  and  cloud-free  regions. 
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Figure  25.  MODIS-derived  cloud-top  heights  - 1825-1830  UTC,  27  Jan  2014. 
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By  combining  the  AIRS  and  MODIS  data  sets,  the  excellent  vertical  resolution  of 
the  AIRS  data  can  refine  the  cloud-top  height  assignment  in  the  MODIS  data  while 
retaining  the  MODIS  1-km  horizontal  resolution.  The  improved  cloud-top  product  could 
be  used  to  determine  expected  cloud-free  line  of  sight  for  an  operational  high  energy  laser 
(HEL)  weapon  system  at  a  specified  location  and  altitude. 

3.4.3  MODIS  Radiance  Data  as  Verification  Tool  for  Laser  Environmental 
Effects  Definition  and  Reference  (LEEDR)  Development 

An  additional  application  for  the  MODIS  data  processing  techniques  developed  for 
cloud-top  height  determination  is  to  obtain  calibrated  radiance  values  for  each  of  the  36 
MODIS  channels.  This  ability  is  useful  for  ongoing  efforts  to  improve  AFIT’s  LEEDR 
modeling  and  simulation  code.  The  LEEDR  software  package  provides  detailed 
atmospheric  characterization,  allowing  users  the  ability  to  customize  key  parameters  such 
as  humidity  profile  and  aerosol  concentrations.  The  user  is  also  able  to  define  cloud  layer 
heights  and  thicknesses,  and  select  cloud  types.  The  location,  date,  and  time  are  input  to 
determine  sun  angle.  Once  the  atmospheric  conditions  are  entered,  LEEDR  calculates 
wavelength  dependent  transmission,  absorption,  and  scattering  along  any  specified  path 
through  the  atmosphere.  Current  improvements  include  a  more  robust  code  to  model  the 
complicated  process  of  multiple  scattering  within  cloud  layers.  The  ability  to  compare 
LEEDR  modeled  upwelling  radiance  with  calibrated  MODIS  measurements  (available  for 
any  location,  time  of  year,  and  cloud  type)  is  proving  invaluable  during  model  refinement. 
Figure  26  and  Figure  27  show  a  granule  of  MODIS  radiance  measurements,  with  seven 
locations  selected  to  give  a  cross-section  of  cloud  type,  height,  and  layer  thickness. 
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Figure  26.  MODIS  1  km  data  granule.  The  seven  squares  in  the  northwest  corner 
of  the  image  correspond  to  the  seven  spectral  radiance  measurements  plotted  in 
Figure  27,  providing  a  sample  of  several  cloud -top  heights  and  cloud  types.  Data 
collected  15  Dec  2014, 1815-1820  UTC. 


Figure  27.  Calibrated  MODIS  radiance  measurements.  Each  spectral  radiance  line 
corresponds  to  one  of  the  MODIS  fields-of-view  identified  in  Figure  26. 
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3.5  Gap  Filling  with  NWP  models 

The  instruments  onboard  the  polar-orbit  A-Train  satellites  collect  data  for  most 
locations  twice  a  day,  so  there  are  11-13  hour  periods  of  time  with  no  updates.  There  are 
two  ways  to  overcome  this  problem.  One  way  would  be  to  exploit  data  collected  from 
similar  instruments  on  satellites  in  different  orbits.  For  example,  Aqua  is  not  the  only 
satellite  with  a  MODIS  instrument.  The  satellite  Terra  has  an  identical  MODIS  unit,  and 
its  orbit  crosses  the  equator  about  three  hours  earlier  than  Aqua,  giving  many  locations  four 
overpasses  a  day  and  cutting  in  half  the  mean  time  between  data  collection. 

A  better  and  more  operationally  useful  way  to  fill  in  these  temporal  gaps  in  satellite 
coverage  is  to  take  advantage  of  the  output  from  NWP  models.  This  research  focuses  on 
incorporation  of  output  from  the  two  models  described  in  Section  2.4.  The  techniques 
developed  in  this  work  to  calculate  Cn2  values  are  applied  to  temperature,  pressure,  height, 
and  wind  data  obtained  from  the  GFS  model  (both  the  0.5  degree  and  0.25  degree 
horizontal  resolution)  and  the  ECMWF  model.  Both  of  these  global  models  have  the  ability 
to  provide  atmospheric  data  for  any  location  of  interest,  but  there  is  a  big  difference  in  how 
accessible  the  data  are.  GFS  output  is  freely  available  to  the  public  for  download  from 
NCEP.  Access  to  ECMWF  forecasts  and  recent  model  output  requires  an  expensive 
license,  limiting  the  number  of  organizations  that  can  work  with  the  data.  The  approach 
taken  in  this  research  is  to  utilize  the  publically  available  ERA-20C  data  which  allows 
analysis  using  data  with  the  same  coverage  and  resolution  as  the  forecast  data.  The  ERA- 
20C  data  has  the  limitation  that  it  only  covers  years  2010  and  earlier. 

The  primary  advantage  to  working  with  output  from  both  the  GFS  and  ECMWF 
models  is  that  it  allows  comparison  between  a  hydrostatic  and  a  non-hydrostatic  model. 
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Another  reason  is  that  it  allows  investigation  into  the  effect  increased  resolution  has  on 
derived  products.  The  ECMWF  grid  spacing  of  0.125  degrees  is  half  and  one-quarter  of 
the  spacing  between  grid  points  in  the  two  GFS  models.  Vertical  coverage  is  another 
advantage  of  the  ECMWF,  with  the  top  pressure  level  of  0.01  hPa  this  model  generates 
output  from  the  surface  to  heights  above  75  km.  In  this  research,  NWP  data  are  used  in 
two  different  ways.  The  first  method  is  to  use  NWP  data  only  to  fill  in  the  temporal  gaps 
between  satellite  overpass.  In  this  method,  the  C„2  derivation  technique  being  analyzed  is 
applied  to  the  satellite  sounder  data  at  the  timesteps  for  which  it  is  available,  and  to  the 
NWP  data  at  3-hour  intervals  in  between.  This  method  allows  the  creation  of  Cn2  profiles 
that  extend  for  several  days  or  weeks  with  enough  data  points  to  make  useful  comparison 
to  scintillometer  measurements. 

The  second  way  NWP  data  are  used  in  this  research  is  to  apply  the  same  Ct2  and 
C2  derivation  techniques  to  data  obtained  by  satellite  sounder,  NWP  model,  and 
rawinsonde  when  multiple  sources  are  available  for  the  same  location  at  nearly  the  same 
time.  This  allows  analysis  of  the  performance  of  the  C2  derivation  technique  when  applied 
to  data  sets  with  different  spatial  and  temporal  resolution,  and  it  also  allows  for 
determination  of  the  suitability  of  the  various  satellite  sounder  and  NWP  products  to  this 
type  of  optical  turbulence  research. 
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4  Results  and  Analysis 


4.1  Validation  of  Satellite-Measured  Temperatures 

The  technique  for  deriving  winds  through  the  thermal  wind  relationship  as 
described  in  Section  3.2  relies  on  often  very  small  horizontal  temperature  gradients, 
making  it  sensitive  to  any  error  in  the  underlying  temperature  data.  Also,  the  full  wind 
profile  is  constructed  by  adding  each  layer’s  wind  gradients  to  the  previous  layer,  which 
causes  even  small  errors  in  measured  temperature  to  compound  and  at  times  lead  to  large 
errors  in  wind  speed  and  direction.  For  these  reasons,  additional  analysis  of  the  error  in 
the  AIRS  temperatures  is  included  as  part  of  this  research. 

The  error  in  AIRS/AMSU-measured  temperatures  varies  with  height,  amount  of 
cloud  cover,  and  geographic  location  (which  affects  the  IR  and  microwave  sounding 
background).  AIRS  reported  error  is  largest  in  the  upper  atmosphere,  where  the  air  density 
is  much  lower  [56],  but  rawinsonde  observation  (RAOB)  temperatures  are  not  available 
for  comparison  at  these  altitudes,  as  the  weather  balloon  does  not  typically  survive  past  an 
altitude  of  25-30  km.  The  height  dependence  in  the  AIRS-reported  error  can  be  seen  in 
Figure  28,  which  plots  an  AIRS  temperature  profile  against  the  rawinsonde  measured 
temperature  profile  for  the  1200  UTC  release  at  the  Anchorage,  AK  station.  The  AIRS 
profile  is  for  the  field-of-view  closest  to  the  rawinsonde  release  coordinates  as  the  satellite 
passed  over  Anchorage  at  1408  UTC. 


49 


AIRS  Temperature  Profile  AIRS-Reported  Temperature  Error  vs  Height 


Figure  28.  Anchorage,  AK  comparison  of  rawinsonde  observation  and  AIRS 
temperatures-15  Oct  2014, 1408  UTC.  In  the  left  plot,  the  size  of  the  error  bar  at 
each  AIRS  height  indicates  the  magnitude  of  the  temperature  error  reported  with 
the  AIRS  level  2  data.  The  green  dotted  line  shows  the  temperature  measurement 
made  by  a  rawinsonde  released  at  the  Anchorage  International  Airport  two  hours 
after  this  AIRS  overpass.  The  right  plot  shows  the  magnitude  of  the  AIRS -reported 
error  as  a  function  of  height. 

For  the  temperature  accuracy  study,  four  locations  in  Alaska  and  Hawaii  were 
selected  due  to  their  scheduled  0000  UTC  and  1200  UTC  rawinsonde  times  closely 
matching  the  AIRS  overpass  times.  For  each  of  these  locations,  all  AIRS  data  for  the  12 
year  period  2003-2014  were  processed,  and  when  the  overpass  time  for  the  sounding 
location  nearest  to  the  rawinsonde  site  was  within  1  hour  of  the  rawinsonde  release,  the 
temperature  profiles  were  compared.  The  results  of  this  analysis  are  presented  in  Figure 
29  and  Figure  30. 
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Figure  29.  Root-Mean-Square  (RMS)  error  between  AIRS  and  rawinsonde 
measured  temperatures  for  four  locations  selected  for  frequent  AIRS  overpass  near 
time  of  0000  UTC  and  1200  UTC  rawinsonde  release.  For  12  years  of  AIRS  data 
(2003-2014),  each  instance  when  the  AIRS  overpass  time  was  within  1  hour  of  the 
rawinsonde  release  is  included.  This  figure  includes  data  for  7519  match-ups  for 
Lihue,  HI,  7358  for  Hilo,  HI,  6792  for  Anchorage,  AK,  and  6414  for  Fairbanks,  AK. 

The  differences  between  the  AIRS -retrieved  temperatures  and  the  rawinsonde 
measurements  are  greatest  in  the  lower  atmosphere  (below  700  hPa)  and  in  the  upper 
atmosphere  (above  20  hPa).  In  the  region  between  20  hPa  and  700  hPa  the  RMS  error  is 
between  1.0  and  2.5  K  for  each  location. 
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In  the  lower  atmosphere  the  satellite-retrieved  temperatures  were  most  accurate  for 
locations  with  low  or  no  cloud  cover,  as  the  resolution  of  the  AIRS  sounder  outperforms 
the  AMSU  instrument.  A  different  effect  was  observed  in  the  upper  atmosphere  for  the 
two  locations  in  Alaska.  For  heights  above  the  100  hPa  level,  the  temperatures  for  cloud- 
covered  locations  were  a  better  match  to  the  rawinsonde  measurements  than  the  cloud-free 
locations  were.  This  effect  could  be  due  to  low  cloud  cover  providing  a  more  uniform 
background  for  the  IR  sounder  or  it  might  identify  the  presence  of  optically-thin  cirrus 
clouds  which  are  not  recognized  by  the  AIRS  cloud  detection  algorithm.  This  effect  is  not 
present  in  the  Hawaii  locations,  as  much  of  the  AIRS  field-of-view  is  water  when  the 
location  is  cloud-free.  Further  investigation  including  additional  locations  and  cloud 
conditions  is  required  to  investigate  the  cause  of  this  effect. 

Figure  30  identifies  a  different  sounding  background  issue  with  the  AIRS 
temperature  profiles.  In  both  Alaska  locations,  there  is  a  strong  cold  temperature  bias  in 
the  data  for  the  few  pressure  levels  closest  to  the  surface.  In  both  cases  there  are  multiple 
levels  near  the  bottom  of  the  profile  for  which  more  than  half  of  the  sampled  temperatures 
are  colder  than  the  rawinsonde  observation  by  more  than  2  K.  It  seems  this  difference  is 
probably  a  winter  time  effect,  when  at  times  a  very  strong  temperature  inversion  exists 
within  500  m  of  the  surface.  A  2010  paper  [57]  on  these  surface-based  inversions  found 
that  in  Fairbanks  in  December  and  January,  the  air  temperature  at  the  surface  is  on  average 
10  K  colder  than  the  air  temperature  at  a  height  of  1000  m.  Due  to  the  vertical  spread  of 
the  weighting  functions,  the  sounding  channels  near  the  bottom  of  the  vertical  profile 
would  receive  a  contribution  (in  this  case  a  reduction  from  the  typical  amount  of  radiance) 
from  the  cold  surface  and  the  first  few  hundred  meters,  which  causes  problems  for  the 
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process  of  assigning  temperatures  to  the  profile.  This  allows  the  cold  surface  to  drag  down 


the  calculated  temperatures  for  the  neighboring  levels.  If  that  is  the  case,  this  effect  should 


show  up  more  in  winter  than  in  summer,  so  analyzing  the  seasons  separately  would  help 


confirm  or  refute  this  explanation. 
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Figure  30.  Temperature  bias  between  AIRS  and  rawinsonde  measured 
temperatures  for  four  locations  selected  for  frequent  AIRS  overpass  near  time  of 
0000  UTC  and  1200  UTC  rawinsonde  release.  These  comparisons  include  the  same 
2003-2014  AIRS  and  rawinsonde  data  as  Figure  29. 
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4.2  AIRS-Derived  Winds 


4.2.1  Accuracy  of  the  AIRS-Derived  Winds 

A  useful  source  of  measured  vertical  profiles  of  wind  speed  and  direction  for  use 
in  analysis  of  AIRS -derived  wind  profiles  is  the  network  of  rawinsonde  observations  [37]. 
These  weather  balloon  soundings  provide  accurate  wind  data  for  hundreds  of  locations 
each  day  at  0000  and  1200  UTC.  There  are  only  two  bands  of  longitude  values  for  which 
the  twice  daily  AIRS  passes  (early  afternoon  and  early  morning  local  time)  are  close  to  the 
time  of  the  RAOBs.  When  they  coincide,  direct  comparison  of  the  wind  profiles  can  be 
made,  as  shown  in  Figure  31. 
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Figure  31.  Comparison  of  AIRS-derived,  RAOB-measured,  and  GFS-modeled  wind 
profiles,  Lihue,  HI,  5  August,  2014,  0000  UTC. 
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Figure  3 1  also  includes  the  corresponding  NWP  wind  profile,  in  this  case  from  the 
Global  Forecast  System  (GFS)  0.5  degree  model.  This  model  is  used  because  both  the 
archived  data  and  forecast  data  are  widely  available,  making  it  a  good  choice  for  research 
and  potential  operational  applications.  GFS  wind  profiles  do  not  match  the  vertical 
resolution  of  the  RAOB,  but  typically  provide  an  accurate  representation  of  the  trends  in 
wind  speed  and  direction.  This  makes  it  practical  to  use  the  GFS  winds  to  evaluate  the 
accuracy  of  the  AIRS-derived  winds,  investigating  locations  and  times  of  day  for  which  no 
RAOB  data  are  available. 

A  grid  of  210  locations  was  selected  in  the  Northern  Hemisphere,  covering  the 
eastern  Pacific  Ocean  and  western  portion  of  North  America  (at  5  degree  latitude  and 
longitude  intervals).  The  geographic  coverage  of  these  sample  sounding  locations  is  shown 
in  Figure  32. 


Coordinates  Selected  for  Evaluation  of  AIRS  Windfinding  Algorithm 


Figure  32.  AIRS  sounding  coordinates  selected  for  GFS  comparison. 
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For  each  of  the  locations  identified  in  Figure  32, 100  AIRS  soundings  were  selected 


for  download,  one  day  and  one  night  sounding  per  week  throughout  2014.  Corresponding 
GFS  data  for  the  nearest  GFS  grid  point  and  closest  3  hour  interval  were  retrieved  for 
comparison.  Of  the  possible  21,000  matches,  18,747  had  both  AIRS  and  GFS  profiles 
available.  These  match-ups  were  used  to  investigate  the  conditions  under  which  the  AIRS 
windfinding  technique  worked  best  and  which  conditions  led  to  the  least  accurate  winds. 

When  calculating  a  wind  profile,  if  the  coordinates  of  interest  fall  near  the  edge  of 
the  AIRS  swath,  the  temperature  and  pressure  gradients  used  in  the  wind  derivation 
technique  may  be  less  accurate  because  of  missing  neighboring  AIRS  data  points.  For  this 
reason,  it  is  expected  that  the  AIRS  winds  should  not  match  the  GFS  winds  as  well  when 
the  closest  AIRS  field-of-view  is  near  an  edge.  Figure  33  shows  the  average  of  all  wind 
speed  profile  RMS  errors,  plotted  by  corresponding  AIRS  field-of-view  number. 
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Figure  33.  Wind  speed  profile  error  by  AIRS  field-of-view.  The  AIRS  scan 
geometry  is  a  left  to  right,  whisk  broom  type,  with  each  scanline  comprised  of  30 
fields-of-view,  numbered  1  through  30. 
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Field-of-view  numbers  1  and  30  represent  the  left  and  right  edges  of  the  AIRS 
swath.  AIRS  wind  speed  was  used  for  this  comparison  because  it  is  more  important  for 
turbulence  derivation  than  wind  direction,  but  similar  results  are  observed  for  u  and  v  wind 
components  individually.  The  least  accurate  results  occur  at  the  edges,  but  it  appears  that 
only  the  outer  two  fields-of-view  are  significantly  impacted  by  this  effect. 

Time  of  day  and  time  of  year  are  additional  factors  that  were  investigated  for  their 
effect  on  the  accuracy  of  the  derived  wind  profile.  Figure  34  and  Figure  35  show  the 
average  RMS  error  values  for  all  of  the  retrieved  soundings,  separated  by  month  and 
whether  it  was  the  local  afternoon  sounding  (day)  or  just  after  local  midnight  sounding 
(night).  The  over  land  soundings  do  not  indicate  any  clear  dependence  on  time  of  day  or 

time  of  year. 
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Figure  34.  Average  AIRS-derived  wind  speed  RMS  error  for  selected  over-land 
soundings,  plotted  by  month. 
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The  over  water  comparisons  in  Figure  35  appear  to  reveal  some  seasonal  variation 


in  accuracy,  evident  in  both  the  day  and  night  sounding  data.  The  RMS  wind  speed  errors 
are  greatest  around  April/May  and  smallest  during  August/September. 
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Figure  35.  Average  AIRS-derived  wind  speed  RMS  error  for  selected  over-water 
soundings,  plotted  by  month. 

The  most  significant  factor  in  how  well  the  AIRS-derived  wind  profile  agrees  with 
the  GFS  wind  is  cloud  cover.  The  atmospheric  temperatures  below  clouds  are  not  retrieved 
by  the  IR  AIRS  instrument — they  are  retrieved  at  reduced  vertical  resolution  by  the  AMSU 
radiometer.  These  AMSU  temperatures  are  used  to  complete  the  temperature  profiles,  and 
the  reduced  vertical  resolution  reduces  the  accuracy  of  the  satellite-derived  winds. 

The  AIRS  reported  cloud  fraction  for  the  soundings  sampled  ranged  from 
completely  clear  (8.1%)  to  completely  overcast  (5.9%).  To  show  this  range,  and  the 
distribution  of  these  cloud  fraction  values,  a  histogram  is  included  as  Figure  36. 
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AIRS  Reported  Field-of-View  Cloud  Fraction 

Figure  36.  Histogram  of  AIRS  reported  cloud  fraction  values,  the  total  effective 
cloud  cover  for  all  cloud  layers  within  the  AIRS  field -of- view. 

Each  AIRS  field-of-view  is  comprised  of  9  separate  IR  soundings,  and  the  total 
effective  cloud  fraction  over  all  cloud  layers  for  these  9  spots  is  determined  by  multiplying 
the  cloud  fraction  for  each  field-of-view  by  the  cloud  emissivity.  The  assumption  is  made 
that  the  clouds  are  spectrally  flat,  and  for  a  cloud  that  is  partially  transmissive,  the 
equivalent  opaque  fraction  is  reported.  This  allows  an  effective  cloud  fraction  value  for 
each  AIRS  field-of-view  to  have  any  value  between  0  and  1.0  [54]  [58]  [59]. 

As  described  in  Section  4.1,  cloud  cover  increases  the  error  in  the  retrieved 
temperatures  in  the  lower  atmosphere.  The  increased  error  in  the  temperatures  retrieved 
near  the  surface  affects  the  accuracy  of  the  winds  in  the  lower  atmosphere,  which  (because 
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of  the  way  the  wind  profiles  are  constructed  from  the  bottom  up)  increases  the  error 
throughout  the  profile. 

This  research  shows  that  the  portion  of  the  AIRS  field-of-view  that  is  cloud 
obscured  affects  wind  profile  accuracy,  with  the  worst  average  accuracy  occurring  when 
cloud  fraction  exceeds  a  value  of  0.8.  Figure  37  shows  the  average  of  all  wind  speed  profile 
RMS  errors  within  each  bin  plotted  against  corresponding  cloud  cover  percentage,  in  5% 
increments. 
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Figure  37.  Average  wind  speed  profile  RMS  error  by  AIRS  reported  cloud  cover 
percentage  for  the  18,747  sampled  AIRS  soundings  and  corresponding  GFS  data 
points. 

Figure  38  depicts  the  spread  of  these  RMS  errors  throughout  the  range  of  cloud 
cover  values,  showing  that  the  cloud  cover  effect  on  wind  profile  accuracy  is  small 
compared  to  the  variability  from  one  sample  to  the  next. 
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Figure  38.  Wind  speed  profile  RMS  error  by  AIRS  reported  cloud  cover  percentage 
for  the  18,747  AIRS  soundings  sampled. 

Another  important  factor  affecting  the  accuracy  of  the  AIRS -derived  wind  profile 
is  the  surface  type  below  the  sounding.  The  AIRS  weighting  functions  for  many  of  the 
lowest  altitude  channels  intersect  the  surface,  and  the  emissive  properties  of  the  surface 
determine  the  amount  of  radiance  contributed  by  the  surface  to  that  channel.  Cloud  cover 
has  a  larger  effect  on  the  AIRS/AMSU  data  over  land  than  over  water,  because  the  included 
AMSU-derived  data  is  more  strongly  affected  by  the  varying  surface  emissivities  found 
over  land.  Figure  39  shows  that  even  when  corrected  for  variations  in  cloud  cover,  the 
surface  emissivity  variation  leads  to  reduced  accuracy  in  regions  where  the  sounder 
background  is  land. 
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Figure  39.  Average  AIRS-derived  wind  speed  profile  RMS  error  (in  m/s),  compared 
with  GFS  wind  profile  for  nearest  grid  point  and  time.  In  this  figure,  the  RMS  error 
for  the  wind  profiles  for  each  of  the  5152  cloud-free  or  nearly  cloud-free  (AIRS 
cloud  fraction  across  the  field-of-view  less  than  0.2)  AIRS  soundings  within  each  5 
degree  by  5  degree  box. 

In  Figure  39,  only  soundings  with  less  than  0.2  cloud  fraction  are  included.  For 
these  soundings,  the  average  wind  speed  RMS  error  for  the  ocean  only  regions  is  8.6  m/s. 
When  land  fills  part  of  the  box,  the  average  RMS  error  is  11.2  m/s,  30%  larger. 

This  technique  of  deriving  winds  from  satellite  sounding  data  has  its  operational 
advantages  (passive  sensing,  single  data  source,  and  computationally  inexpensive),  but  the 
accuracy  of  the  resulting  winds  does  not  match  the  winds  available  from  the  NWP  models 
investigated.  When  compared  with  rawinsonde  measured  winds,  the  AIRS  wind  profiles 
obtained  often  capture  trends  in  wind  speed  and  direction,  but  the  satellite-derived  wind 
profiles  are  typically  not  as  accurate  as  winds  provided  by  modern  NWP  models. 
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4.2.2  AIRS-Derived  Winds  in  the  Stratosphere 

While  not  found  to  be  suitable  for  low-altitude  operational  applications,  the  AIRS- 
based  windfinding  technique  can  provide  wind  information  in  regions  of  the  atmosphere 
not  currently  covered  by  the  GFS  model.  The  AIRS  data  retrieves  temperatures  up  to 
altitudes  of  70-75  km  (standard  pressure  level  0.0161  hPa),  so  derived  wind  profiles  extend 
throughout  the  stratosphere  and  lower  mesosphere,  above  the  highest  data  available  from 
most  current  NWP  models.  Direct  measurements  of  winds  in  the  upper  atmosphere 
(heights  unreachable  by  weather  balloons)  are  infrequent  and  expensive.  The  ability  to 
diagnose  winds  at  these  altitudes  from  satellite  sounder  data  would  be  operationally  useful. 

The  traditional  way  to  obtain  winds  at  altitudes  between  40  and  80  km  is  to  launch 
a  rocketsonde.  An  example  of  this  is  illustrated  in  Figure  40,  which  shows  the  full  AIRS- 
derived  u  and  v  wind  profiles  compared  with  wind  components  measured  by  a  series  of 
four  Super-Loki/PWN  12A  rockets  deploying  Rocket  Balloon  Instrument  (ROBIN) 
Inflatable  Falling  Spheres.  At  altitude,  each  ROBIN  sphere  is  inflated  by  evaporation  of  a 
contained  liquid.  It  has  an  attached  corner  reflector,  allowing  it  to  be  tracked  by  a  precision 
ground-based  radar  as  it  falls  through  the  atmosphere  [60].  While  it  does  not  match  the 
vertical  resolution  of  the  rocketsonde  data  in  the  stratosphere,  the  AIRS  windfinding 
technique  captures  the  wind  component  trends,  providing  some  insight  into  the  strength 
and  direction  of  the  wind  at  heights  well  above  those  modeled  by  GFS  and  other  easily 
accessed  NWP  models. 
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Figure  40.  Comparison  of  the  u  and  v  wind  component  profiles  derived  from  AIRS 
temperature  data  with  winds  measured  by  rawinsonde  and  by  rocketsondes 
launched  from  the  Pacific  Missile  Range  Facility  (PMRF)  on  8  Jun  2015.  The 
rocketsonde  data  was  obtained  through  personal  communications  with  Mary 
Bedrick,  AFLCMC/XZIG,  Wright-Patterson  AFB  and  Dr.  Clara  O’Farrell,  NASA 
JPL  [61]. 

4.3  Calculation  of  Refractive  Index  Structure  Function  ( Cn 2) 

4.3.1  AIRS-Derived  C,r  Profiles 

One  of  the  challenges  with  the  AIRS-based  Cn2  derivation  technique  is  the  result  of 
the  nature  of  the  AIRS  temperature  data.  Temperature  retrieval  by  IR  sounding  involves 
assigning  a  temperature  to  each  altitude  based  on  the  radiance  received  on  many  channels. 
The  weighting  function  for  each  channel  defines  how  sensitive  that  channel  is  to  radiance 
from  each  layer  of  the  atmosphere.  These  weighting  functions  overlap  each  other 
vertically,  and  the  radiation  received  by  an  individual  channel  emanates  from  a  range  of 
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altitudes.  Through  an  iterative  process,  the  temperature  assigned  to  each  level  by  the 
retrieval  algorithm  is  based  on  the  radiance  received  by  many  channels  with  weighting 
function  peaks  near  the  height  of  that  level.  This  leads  to  an  unavoidable  smoothing  of  the 
temperature  data — even  when  the  AIRS  temperatures  match  RAOB  temperatures  within  a 
few  K  at  every  altitude,  the  weather  balloon  is  able  to  detect  temperature  structure  and 
variability  over  small  vertical  distances  that  are  smoothed  out  by  the  satellite  retrieval 
process. 

This  temperature  smoothing  leads  to  less  variability  than  is  present  in  directly 
measured  wind  profiles.  Figure  41  shows  an  example  of  this  smoothing  in  the  left  plot.  In 
this  example,  the  AIRS -derived  wind  speed  agrees  very  well  with  the  rawinsonde 
measurement,  but  the  rawinsonde  measures  speed  fluctuations  of  5-10  m/s  throughout  the 
balloon’s  ascent.  These  fluctuations  are  not  present  in  the  AIRS  wind  profile. 

When  the  AIRS -derived  wind  gradient  is  used  to  calculate  the  gradient  Richardson 
number,  the  smoothing  of  the  wind  speed  profile  leads  to  a  larger  Richardson  number.  The 
right  plot  in  Figure  41  shows  a  comparison  between  the  Richardson  numbers  calculated 
from  the  AIRS  and  the  RAOB  data.  The  AIRS-derived  values  are  between  1  to  2  orders 
of  magnitude  larger  than  the  RAOB  derived  values  for  much  of  the  profile.  This  higher 
gradient  Richardson  number  leads  to  lower  Cn2  values,  especially  at  heights  where  the 
smoothed  wind  speed  reaches  relative  minima  or  maxima.  If  the  wind  speed  does  not 
change  much  from  one  level  to  the  next  (e.g.,  the  levels  between  11  and  12  km  in  Figure 
41)  the  denominator  of  the  Richardson  number  equation  approaches  zero,  leading  to  a  spike 
in  Ri,  and  a  corresponding  proportional  drop  in  Kh/Km  and  C,,2.  This  demonstrates  that 
even  though  the  sounder  data  has  comparable  vertical  resolution  and  good  temperature 
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accuracy  (most  of  the  temperature  profile  within  1.5  K  of  RAOB  measured  temperature), 
there  are  problems  introduced  by  the  fact  that  space-based  sounding  will  always  fail  to 
measure  temperature  structure  and  fluctuations  that  occur  on  a  scale  smaller  than  2  or  3 
vertical  data  points.  Missing  the  temperature  variability  has  an  impact  on  the  wind  profile, 
which  affects  the  Richardson  number  almost  exclusively  in  one  direction — pushing  Ri  to 
larger  values  which  leads  to  smaller  Cn2  values. 


Figure  41.  Smoothing  of  satellite  data  leads  to  large  gradient  Richardson  numbers  - 
11  Apr  2014  1250  UTC,  Anchorage,  AK. 

4.3.2  Use  of  NWP  Winds  with  AIRS  Temperatures  to  Calculate  Cn2 

An  alternative  source  of  wind  information,  which  avoids  the  accuracy  problems 
with  the  AIRS-only  winds,  is  NWP  data.  NWP  winds  match  rawinsonde  measurements 
more  closely  than  the  AIRS-derived  winds  do,  so  their  use  for  C,,2  calculation  is  expected 
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to  yield  better  results.  While  not  capturing  all  of  the  variability  in  component  wind  speeds 
over  small  vertical  distances  that  are  measured  by  a  rawinsonde,  GFS  winds  are  not  as 
smoothed  as  the  AIRS-derived  winds  (see  Figure  3 1  for  an  example).  This  difference  leads 
to  more  structure  in  the  Cn 2  profile,  as  seen  in  Figure  42,  where  the  C,,2  values  calculated 
using  GFS  0.5  degree  model  winds,  obtained  through  the  NOMADS  server,  are  compared 
to  an  AIRS-only  Cn2  profile. 


Figure  42.  Comparison  of  Cn  profiles  from  AIRS-only  technique  and  hybrid  AIRS 
temperatures  and  GFS  winds.  This  comparison  is  for  Lihue,  HI,  2  January  2013, 
1800  UTC. 


67 


4.3.3  Creation  of  a  Standard  Maritime  Cn2  Model 

The  techniques  developed  in  this  research  project  provide  the  tools  necessary  to 
create  a  standard  Cn2  model  for  specific  locations  or  conditions  by  calculating  hundreds  of 
profiles,  across  a  range  of  different  weather  conditions,  and  averaging  them.  For  each 
height,  a  mean  C2  value  and  the  standard  deviation  can  be  determined.  Even  though  C2 
values  at  a  given  location  can  vary  greatly  as  atmospheric  conditions  change,  at  times  by 
more  than  an  order  of  magnitude  from  one  minute  to  the  next,  a  standard  model  is  a  useful 
engineering  tool,  providing  a  way  to  determine  an  envelope  of  expected  performance  for  a 
directed  energy  system  or  sensor  in  a  specified  environment.  The  widely  used  C2  model, 
Hufnagel- Valley  5/7  (previously  presented  in  Figure  11),  was  developed  to  characterize 
optical  turbulence  (C„2)  as  it  varies  with  height  over  land.  While  the  H-V  5/7  model  is  not 
meant  to  predict  precise  values,  observations  will  likely  be  in  a  range  close  enough  to  the 
model’s  prediction  to  make  it  a  useful  system  design  and  planning  tool.  H-V  5/7  is  not 
designed  for,  nor  well-suited  to,  modeling  turbulence  over  the  ocean.  Thus  recent  AFIT 
research,  conducted  by  Ft  Greg  Anderson,  has  focused  on  the  development  of  an  H-V-like 
standard  model  of  optical  turbulence  for  a  maritime  atmosphere  [62],  Maritime 
temperature  climatology  data  collected  by  AIRS  and  wind  data  from  the  GFS  model  were 
used  in  the  Tatarskii  relationship  as  described  in  Section  3.3.1  to  calculate  over  29,000  C„ 2 
profiles  over  the  world’s  oceans.  Three  standard  models  were  proposed  as  alternatives  to 
the  H-V  model.  Results  showed  that  maritime  profiles  generally  do  not  exhibit  the  surface 
spike  in  turbulence  seen  in  H-V  5/7.  Additionally,  a  strong  latitudinal  variation  in  the 
height  of  the  Cn2  inflection  associated  with  the  height  of  the  tropopause  is  observed, 
motivating  the  need  for  separate  models  for  polar  and  tropical  regions. 
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For  this  research  project,  40  maritime  locations  were  selected,  and  for  each  of  those 
locations,  all  available  day  and  night  AIRS  soundings  were  collected  for  the  entire  year  of 
2013.  The  locations  were  selected  to  sample  a  range  of  global  maritime  atmospheres,  and 
allow  day-night,  seasonal,  and  regional  comparisons  to  be  made.  The  40  locations  used  to 
generate  these  standard  Cn 2  profiles  are  shown  in  Figure  43. 


Figure  43.  Satellite  and  NWP  data  were  collected  for  40  global,  maritime  locations, 
identified  by  the  numbered  markers  in  this  figure  [62]. 

A  plot  of  the  calculated  and  median  daytime  C„2  values  for  the  month  of  April  2013, 
for  a  maritime  location  north  of  the  state  of  Alaska  is  shown  in  Figure  44.  The  observed 
values  represent  every  Cn2  value  calculated  from  the  combined  AIRS  and  GFS  wind  data 
available  for  this  location  during  this  time  period. 
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Figure  44.  One  month  of  daytime  AIRS  and  NWP  data  derived  Cn 2,  April  2013.  The 
location  is  identified  as  position  2  in  Figure  43  [62]. 

This  averaging  process  was  expanded  to  include  all  40  maritime  locations  and 
extended  for  the  entire  year  of  2013.  One  example  of  the  resulting  standard  Cn 2  models  is 
a  dynamic  piecewise  model,  which  is  comprised  of  two  Cn 2  profiles  to  account  for  the 
variability  in  tropopause  height  with  latitude.  Figure  45  shows  this  model  (for  polar  and 
tropical  regions),  each  plotted  over  the  AIRS  and  NWP  derived  Cn2  values  for  all  of  the 
selected  locations  and  the  entire  2013  time  period  of  the  study. 


70 


Piecewise  Fit,  Polar  Piecewise  Fit,  Tropic 


Figure  45.  Dynamic  piecewise  polar  and  tropical  model  Cn2  profiles  [62]. 

The  techniques  used  to  create  this  maritime  model  rely  on  meteorological  satellite 
and  NWP  data,  exploiting  the  combination  of  AIRS  atmospheric  temperature  data  and 
NWP  wind  profiles.  This  combination  provides  an  excellent  data  set  for  creating  this  type 
of  standard  Cn 2  profile,  exploiting  the  vertical  resolution  of  the  AIRS  instrument  with  the 
accurate  horizontal  wind  speeds  of  an  NWP  model,  GFS  in  this  case.  By  further  filtering 
the  AIRS  data,  selecting  a  smaller  geographic  window,  narrowing  to  a  specific  time  of 
year,  or  separating  day  and  night  overpasses,  this  technique  could  be  employed  to  develop 
custom,  operationally  relevant  standard  Cn 2  models  to  meet  specific  requirements.  With 
over  twelve  years  of  AIRS  data  available,  it  is  possible  to  obtain  data  from  enough 
overpasses  to  generate  a  statistically  significant  sample  size  for  narrow  conditions.  For 
example,  a  standard  model  generated  using  only  daytime,  winter  soundings  within  100  km 
of  a  single  midlatitude  location  could  still  include  data  from  more  than  80,000  AIRS 
soundings  during  the  2002-2015  time  period. 
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4.3.4  Derived  Cn2— NWP-Only 

The  availability  of  NWP  data  at  three  hour  intervals,  or  even  one  hour  intervals  for 
some  models  recently  put  into  operation,  assures  that  for  operational  requirements  below 
30  km  in  height,  forecast  data  are  available  within  90  minutes  of  any  time  of  interest.  NWP 
models  provide  all  of  the  atmospheric  parameters  required  to  calculate  Cn2  using  the 
traditional  Tatarskii  approach.  One  advantage  of  this  NWP-only  technique  is  the  ability  to 
forecast  Cn2  for  a  specified  location  in  the  future,  something  not  possible  if  the  calculation 
cannot  begin  until  30-60  minutes  after  the  satellite  passes  overhead  and  the  temperature 
data  are  processed. 

The  Cn2  forecast  could  extend  as  far  into  the  future  as  the  underlying  NWP  model 
(10-16  days  for  GFS  and  7  days  for  ECMWF),  but  of  course  the  potential  error  in  the 
forecast  increases  with  the  length  of  the  forecast.  Figure  46  shows  example  vertical  Cn2 
profiles  generated  solely  from  NWP  data.  In  this  example,  ECWMF  0.125  degree 
horizontal  resolution  data  for  a  20  x  30  degree  geographic  region  covering  the  southeast 
United  States  are  used  to  calculate  both  a  profile  for  a  single  grid  point  and  a  mean  profile 
for  all  profiles  within  the  region. 


72 


NWP-Derived  C2,  Single  Profile  NWP-Derived  C2,  Mean  Profile 

n’  a  n’ 


Figure  46.  Cn2  profiles  for  sample  data  set  from  the  ECMWF  ERA-20C  re-analysis 
model  for  1800  UTC,  1  Jun  2010.  The  left  plot  is  for  a  single  grid  point  at  the 
northwest  corner  of  the  sampled  data  (40°  N,  95°  W),  and  the  right  plot  is  the  mean 
Cn  value  at  each  height  for  the  38,400  vertical  profiles  across  the  southeast  United 
States  and  Gulf  of  Mexico  (the  geographic  region  depicted  in  the  top  right  corner  of 
the  plot). 

Figure  47  shows  the  range  of  derived  Cn  values  throughout  this  ECMWF  data  set. 
Throughout  the  vertical  profile,  80%  of  the  Cn  values  fall  within  an  envelope  from  about 
1.5  orders  of  magnitude  below  to  about  1  order  of  magnitude  above  the  median  value. 
These  38,400  profiles  are  all  for  the  same  time,  for  a  20°  band  of  latitude,  so  they  exhibit 
some  common  characteristics.  For  example,  the  bump  in  Cn  present  in  nearly  all  of  the 
profiles  between  13  and  15  km  is  indicative  of  the  height  of  the  tropopause  at  this  time  of 
year  in  this  geographic  region. 
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Figure  47.  Cn2  profiles  for  38,400  grid  points  within  a  southeast  US  data  set  from 
the  ECMWF  ERA-20C  re-analysis  model  for  1800  UTC,  1  Jun  2010.  The  top  plot 
contains  the  entire  vertical  data  set.  The  bottom  plot  is  the  bottom  40  km  of  the 
same  data  set  to  better  show  the  spread  in  values  around  the  median  profile. 
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4.3.5  Modifications  to  Cn2  Calculation  Technique 


4.3.5.1  Additional  Wind  Shear  Terms  in  Richardson  Number 

Turbulent  Kinetic  Energy  (TKE)  is  the  kinetic  energy  associated  with  turbulent 
fluctuations  of  the  atmosphere.  A  useful  quantity  in  the  study  of  turbulence  and  dynamic 
meteorology  is  the  turbulent  kinetic  energy  per  unit  mass  [26] 
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(27) 


The  rate  of  change  of  TKE/m  is  determined  by  all  sources  or  losses  of  turbulence  within 


the  volume  of  interest.  These  generation  and  loss  terms  are  broken  out  in  Equation  (28): 


d (TKE  in) 

— - - — -  =  Ad  +  M  +B  +  Tr-s  (28) 

dt 


where  Ad  is  the  advection  of  turbulence  by  the  mean  wind,  M  is  mechanical  generation, 
and  B  is  buoyant  generation  or  consumption,  Tr  is  transport  by  the  turbulence  itself,  and  s 
is  the  term  for  viscous  dissipation  [26].  Of  these  terms,  the  two  most  important  to  this 
research  are  mechanical  generation  and  buoyant  generation  or  consumption  [63],  One  way 
to  represent  the  balance  of  these  is  the  gradient  Richardson  number  in  wind  component 
form, 


The  denominator  of  Equation  (29)  represents  the  mechanical  generation  of  turbulence  by 
wind  shear,  where  the  two  terms  represent  the  vertical  shear  in  the  two  components  of  the 
horizontal  wind,  u  and  v.  These  terms  are  the  only  two  wind  shear  terms  available  from  a 
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single  vertical  wind  profile,  such  as  the  set  of  wind  measurements  made  by  a  rawinsonde. 
There  are  additional  wind  shear  terms  that  could  be  responsible  for  generation  of 
turbulence.  The  full,  3-dimensional  mechanical  generation  term  could  be  expressed  as 


While  the  additional  components  are  not  available  from  a  single  rawinsonde  measurement, 
the  3-D  volumes  of  wind  data  studied  in  this  research  does  provide  the  ability  to  calculate 
them.  The  wind  field  calculated  from  applying  the  AIRS  windfinding  technique  to  a  series 
of  neighboring  AIRS  fields-of-view  provides  the  horizontal  shear  in  the  horizontal  wind 

terms,  —  and  — .  Vertical  winds  are  not  available  from  the  AIRS  windfinding 

dy  dx 

technique,  but  they  are  available  from  NWP  models. 

3-D  numerical  weather  models  contain  vertical  wind  components  at  each  grid  point 
by  tracking  convergence  and  divergence  of  the  modeled  horizontal  winds.  Starting  from 
the  surface  model  level,  horizontal  winds  converging  at  a  grid  point  force  air  up  (because 
it  cannot  go  into  the  ground)  and  diverging  horizontal  winds  pull  air  down  from  above. 
Once  all  of  the  vertical  motion  is  calculated  for  the  lowest  level,  these  vertical  wind 
components  specify  the  bottom  boundary  condition  for  the  next  higher  level.  Using  the 
divergence  and  convergence  of  the  horizontal  winds  at  the  second  level,  that  level’s  vertical 
motion  can  be  calculated,  and  so  on  throughout  the  3-D  grid. 

As  highlighted  in  Section  2.4.2,  when  compared  to  a  hydrostatic  model,  a  non¬ 
hydrostatic  model  more  accurately  represents  vertical  motion  because  in  addition  to 
vertical  motion  forced  by  convergence  and  divergence,  buoyancy  forces  are  modeled,  and 
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the  vertical  accelerations  calculated  allow  the  winds  to  evolve  in  a  more  accurate  manner 


from  timestep  to  timestep.  An  approach  taken  in  this  research  is  to  exploit  the  wind  fields 
available  from  the  non-hydrostatic  ECMWF  model  to  calculate  a  modified  gradient 
Richardson  number,  using  all  six  of  these  shear  terms,  to  investigate  the  impact  of  the 
inclusion  of  the  additional  components. 

First  for  a  single  sample  data  set,  the  original  and  modified  gradient  Richardson 
number  profiles  are  calculated  and  compared.  The  effect  of  the  additional  components  is 
very  small.  Figure  48  shows  the  difference  between  the  original  and  modified  Richardson 
number  profiles  for  this  test  case. 


Figure  48.  Comparison  between  gradient  Richardson  number,  and  modified  3-D, 
six  shear  term  gradient  Richardson  number,  which  accounts  for  additional  sources 
of  mechanical  generation  of  turbulence.  This  profile  is  calculated  for  a  sample  set  of 
coordinates  40  km  east  of  Springfield,  MO,  on  1  June  2010  at  0000  UTC.  The  two 
largest  percent  decreases  for  this  example  are  16%  and  8%,  and  at  all  other  heights 
the  difference  is  less  than  2%.  The  mean  difference  is  0.0032%. 
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The  slight  decrease  in  Richardson  number  leads  to  an  increase  in  the  ratio  of 
Kh /Km  and  a  corresponding  increase  in  C„2,  but  the  increase  would  be  slight  (much  less 

than  1%)  throughout  most  of  the  vertical  profile.  The  condition  that  would  likely  allow  the 
modification  to  have  a  significant  effect  is  a  height  where  there  is  almost  no  change  in  the 
horizontal  wind  components  from  one  level  to  the  next.  In  Figure  49  the  vertical  profiles 
of  the  six  shear  components  are  plotted  individually  for  the  previous  wind  profile. 
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Figure  49.  Six  wind  shear  components  that  make  up  the  mechanical  generation  of 
turbulence  term  in  the  modified  3-D  gradient  Richardson  number.  These  profiles 
are  calculated  for  the  sample  set  of  coordinates  detailed  in  Figure  48. 

Throughout  much  of  the  profile,  the  two  shear  terms  included  in  the  traditional 
gradient  Richardson  number  dominate,  typically  2-3  orders  of  magnitude  larger  than  the 
components  added  by  the  modification.  Figure  49  only  represents  the  shear  profiles  for  a 
single  vertical  path  through  the  ECMWF  data.  To  better  understand  the  typical  differences 
in  magnitude  between  the  six  components,  Figure  50  reproduces  the  previous  plot,  but  for 
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the  mean  at  each  height  from  many  profiles.  A  160  x  240  grid  was  assembled  from  the 
0.125  degree  horizontal  resolution  ECMWF  ERA-20C  re-analysis  data  set.  For  each  of 
these  38,400  grid  points,  the  six  shear  components  for  the  91  vertical  levels  are  calculated, 
and  the  profiles  of  the  mean  values  of  each  component  are  plotted  in  Figure  50.  These  are 
not  independent  samples,  as  they  are  adjacent  coordinates  in  the  same  data  set,  but  the  area 
is  large  enough  to  cover  a  range  of  different  weather  conditions,  giving  a  good  indication 
of  where  the  mean  magnitudes  would  settle  out  as  additional  data  from  different  dates  and 
times  are  included. 
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Wind  Shear  Components  (s'1) 

Figure  50.  Average  values  for  each  of  the  six  wind  shear  components  that  make  up 
the  mechanical  generation  of  turbulence  term  in  the  modified  3-D  gradient 
Richardson  number.  These  profile  averages  are  calculated  for  a  0.125  degree 
ECMWF  ERA-20C  grid  between  20-40°  N  latitude  and  65-95°  W  longitude,  for  1 
Jun  2010. 
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Comparison  of  the  average  profiles  confirms  the  trend  observed  in  Figure  49’ s 


0  —  q- 

single  profile.  The  two  shear  components  —  and  —  are  the  most  significant.  Next  are 

dz  8z 


the  horizontal  shear  in  the  horizontal  wind  terms  —  and  — ,  and  because  the  vertical 

dy  ax 


wind  speeds  are  so  much  smaller  than  typical  horizontal  wind  speeds,  the  components  of 

dw  dw 

shear  in  the  vertical  wind  —  and  —  are  the  least  significant.  To  further  minimize  the 

dy  dx 


impact  of  the  four  smaller  shear  terms,  these  gradients  are  each  squared  when  included  in 
the  Richardson  number  calculation.  While  there  is  a  physical  basis  for  including  these 
additional  shear  components,  the  accuracy  and  usefulness  of  the  resulting  C„ 2  values  are 
not  significantly  improved  when  compared  to  the  traditional  approach  of  only  considering 
the  two  dominant  terms. 


4.3.5.2  Problem  with  Tatarskii  technique  in  a  well-mixed  layer 

The  traditional  Tatarskii  approach  to  calculating  Cn 2  has  an  inherent  problem  in  a 
well-mixed,  neutrally-buoyant  layer.  Cr  is  proportional  to  the  local  deviation  from  the  dry 
adiabatic  lapse  rate.  In  a  neutrally-buoyant  layer  where  the  actual  lapse  rate  is  very  close 
to  the  dry  adiabatic  lapse  rate,  Ct,  and  therefore  Cn2,  can  be  vanishingly  small.  This  small 
value  can  underrepresent  the  strength  of  the  turbulence  present,  as  the  assumption  that  all 
turbulence  is  related  to  local  buoyant  generation  neglects  additional  sources.  Turbulence 
can  be  generated  by  wind  shears,  leading  to  mixing  which  tends  to  push  the  lapse  rate 
toward  adiabatic  so  the  presence  of  the  turbulent  mixing  would  not  be  reflected  in  the 
vertical  temperature  gradient.  Also,  turbulence  may  be  advected  into  the  area  of  interest 
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by  the  mean  wind,  where  it  may  persist  for  a  time  even  if  the  local  vertical  temperature 
gradient  is  not  favorable  to  buoyant  generation. 


4.3.5.3  Inclusion  of  Pressure  Gradient  Term 


The  traditional  Tatarskii  approach  used  in  this  research  (detailed  in  Section  3.3.1) 
derives  Cn2  from  the  temperature  structure  function  Cr2  and  dn/ dT  according  to 


C2  = 


dn 

ySTj 


c2 


(31) 


which  led  to  Equation  (23).  Tatarskii  [35]  makes  the  assumption  that  the  turbulence  is 
incompressible  and  that  for  optical  wavelengths  in  dry  air,  gradients  in  refractive  index 
depend  only  on  the  vertical  potential  temperature  gradient.  With  this  assumption  removed, 
Cn  can  be  calculated  by  [64] 


c;=* 


f  TS  \ 
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L  % 


dn  dO  dn  dp'  dn  dev 
■  + - — +  -  v 


dT  dz  dp  dz  dev  dz 


(32) 


The  potential  temperature  gradient  term,  —  — ,  is  the  dominant  term  and  does  lead  to  the 

dT  dz 

largest  contribution  to  the  strength  of  the  optical  turbulence,  but  the  additional  terms  are 
non-zero  and  in  some  cases  are  important  to  include.  At  optical  wavelengths,  the  vapor 
dn  de„ 


pressure  term 


de  dz 


is  still  small  enough  to  be  neglected  in  this  research,  but  the 


pressure  term  ^  can  be  an  important  source  of  optical  turbulence,  especially  in  well- 

mixed,  neutrally-buoyant  atmospheres  where  the  potential  temperature  gradient  term 
becomes  very  small.  This  pressure  gradient  term  is  discussed  in  detail  in  Section  4.3.5. 5. 
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Data  sets  generated  by  high  resolution,  non-hydrostatic  NWP  models,  such  as  the  ECMWF 
ERA-20C  atmospheric  re-analysis  model,  allow  comparison  of  the  magnitude  of  each  of 
the  terms  in  Equation  (32)  for  various  atmospheric  conditions. 


4.3.5. 4  Pressure  Fluctuation  Effect  on  Index  of  Refraction 

The  index  of  refraction  of  air  depends  on  the  local  temperature,  pressure,  and 
humidity,  so  the  basic  equation  of  state  is  [22]  [35] 

N  -  N(T ,  p,q)  (33) 

where  N  is  related  to  the  refractive  index  by  N  =  (n-l)xlO6  ,  T  is  temperature,  p  is 
pressure,  and  q  is  absolute  humidity  (mass  of  water/volume).  Equation  (33)  is 
differentiated  to  yield 


7AT  y  1rT1  y  1  y  7 

dN  =  —  dl  +  —  dp  +  —  dq 


To  determine  Cn2,  the  local  fluctuations  in  refractive  index  are  approximated  by 


dN)  0  +  (dN_ 


,  dN 

p  +  — 

(dq 


where  9,  p' ,  and  q  are  the  fluctuations  about  the  mean  value  for  each  quantity. 


As  experimentally  determined  and  published  by  Barrell  and  Sears  in  1939  [65], 
Equation  (36)  provides  a  way  to  calculate  the  index  of  refraction  of  air  as  a  function  of 
temperature,  pressure,  and  humidity. 
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(ft-l)xlQ6  =|  0.378125 +  Q'QQ2|414  +  Q'QQQ(y 794 
_  A  l4 

p[l  +  (l.049-0.157r)/7xl0'6] 
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1  +  77273.16 


0.000680 
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/ 


1  +  77273.16 


(36) 


The  numerical  coefficients  in  this  equation  are  appropriate  for  their  chosen  units, 
where  X  is  wavelength  in  pm,  p  is  pressure  in  mm  Hg,  T  is  temperature  in  degrees  Celsius, 
and/is  the  partial  pressure  of  water  vapor  in  torr  [65],  To  simplify  Equation  (36)  the  two 
coefficients  that  define  the  wavelength  dependent  dispersion  are  redefined  as  A(X)  and  B(X) 
[22]: 

nrtno«~K  0.0021414  0.00001794 

A(A)  =  0.378125  H - 2 - H - 4 -  (37) 


and 


B(y l)  =  0.0624- 


0.000680 

/l2 


(38) 


The  values  of  these  dispersion  coefficients  throughout  the  visible  spectrum  and  into  the 
near-IR  are  plotted  in  Figure  51.  The  effect  of  these  dispersion  coefficients  across  the 
visible  and  near  IR  spectrum  is  small  [22],  and  in  this  research  a  characteristic  wavelength 
near  the  middle  of  the  range  of  wavelengths  of  interest  will  be  selected  and  used  throughout 
the  refractive  index  calculations  that  follow.  The  wavelength  selected  is  0.5  pm,  yielding 
the  values  A(X)  =  0.387  and  B(X)  =  0.0597.  It  is  noted  that  the  error  introduced  by  this 
choice  is  1.35%  at  300  nm,  and  1.7%  at  1.0  pm. 
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Figure  51.  Barrell  and  Sears  dispersion  coefficients. 

With  the  simplification  by  replacing  the  dispersion  coefficients,  Equation  (36) 
becomes 


(n-l)xlO6  =  A(A) 


p[l  +  (l.049-0.157r)pxl0“6] 


1  +  7/273.16 


-BU) 


f 


1  +  7/273.16 


(39) 


Changing  the  numerical  coefficients  to  be  consistent  with  the  units  used  throughout  this 
document,  replacing  (n-l)xlO6  with  N,  and  setting /to  zero  for  dry  air,  Equation  (39)  is 
modified  to  become 


Ndry  =  A(A) 


204.9 p 


1  +  i?(3.295xl0“5  — 1.1 78x1 0-7  7) 


(40) 


where  p  is  now  pressure  in  hPa  and  7  is  temperature  in  K.  It  is  this  form  of  the  equation 
that  is  differentiated  to  determine  the  effect  that  the  temperature  and  pressure  fluctuations 
have  on  refractive  index  according  to  Equation  (35).  Again,  the  effect  of  humidity  is  small 
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at  the  wavelengths  of  interest,  and  the  third  term  in  Equation  (35), 


rdN^ 
ydQjjj.Q 


q  ,  is 


neglected  in  this  research. 

The  partial  derivatives  of  Equation  (40)  with  respect  to  temperature  and  pressure 


are 


and 
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4.3.5.5  Cn  Values  Calculated  using  Temperature  and  Pressure  Gradients 

The  partial  derivatives  are  entered  into  this  dry  air  expression  for  Cn 


Cl  =  a1 
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(43) 


The  term  —  is  the  change  in  refractive  index  with  change  in  temperature  and  —  is  the 
dT  dp 

change  in  refractive  index  with  change  in  pressure,  both  determined  by  differentiating  the 

d0 

Barrell  and  Sears  equation,  as  shown  in  Equations  (41)  and  (42).  The  term  —  is  the 

dz 

deviation  from  dry  adiabatic,  as  previously  calculated  for  the  traditional  Tatarskii  approach. 

The  interesting  term  is  the  local  vertical  pressure  fluctuation  term,  .  Here  p'  is  not 

dz 

atmospheric  pressure,  but  the  fluctuation  in  atmospheric  pressure.  The  approach  taken  in 
this  research  is  to  follow  the  train  of  thought  that  if  the  atmosphere  were  at  all  times,  and 
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at  all  locations,  in  hydrostatic  balance  there  would  be  no  local  vertical  pressure  fluctuations. 
Furthermore,  if  the  pressure  at  one  location  were  perturbed  from  hydrostatic  balance,  it 
would  be  a  transient  event  as  the  surrounding  atmosphere  would  immediately  apply  a 
restoring  force  to  move  back  toward  hydrostatic  balance.  This  is  an  argument  for  stating 
that  the  magnitude  of  any  deviations  from  hydrostatic  equilibrium  are  of  the  same  order  as 
the  transient  pressure  fluctuations.  Therefore,  if  there  were  a  way  to  identify  at  one  instant, 
at  one  location,  the  local  pressure  deviation  from  hydrostatic  equilibrium,  that  deviation 
could  be  used  to  estimate  the  magnitude  of  the  temporal  pressure  fluctuation  at  that 

location,  dp' .  Therefore  is  the  vertical  pressure  gradient  deviation  from  hydrostatic. 

dz 

A  non-hydrostatic  model,  like  the  ECMWF,  makes  this  possible.  The  technique 
begins  with  bottom  level  of  the  model  and  working  up  one  level  at  a  time  from  there,  using 
a  rearranged  hypsometric  equation: 


P  2 


Pi 


g&z/RdTv 


(44) 


At  model  level  2,  this  equation  is  used  to  calculate  what  the  pressure  would  be  if  hydrostatic 
forces  were  in  balance,  using  the  model’s  pressure  at  level  1,  the  difference  between  the 
heights  of  levels  1  and  2,  and  the  mean  virtual  temperature  in  the  layer  between  levels  1 
and  2.  The  difference  between  what  the  pressure  would  be  (if  exactly  in  hydrostatic 
balance)  and  what  the  pressure  at  level  2  actually  is  in  the  model  defines  the  local  pressure 
deviation  dp' .  This  pressure  deviation  is  divided  by  the  vertical  distance  between  the 


7  f 

levels  to  approximate - .  This  process  is  repeated  for  each  higher  grid  point  in  the  model 

dz 
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to  calculate  a  pressure  deviation  value  for  every  level  in  the  model  except  the  surface  level. 
ctjj f 

These  values  for  —  are  used  in  Equation  (43)  to  calculate  C2  values  that  are  based  on 

dz 

the  temperature  deviations  from  adiabatic  and  the  pressure  deviations  from  hydrostatic. 

To  illustrate  the  effect  of  this  modification,  a  vertical  ECMWF  ERA-20C  profile  is 
selected  because  it  contains  a  layer  between  600  and  800  m  where  the  lapse  rate  very  nearly 
matches  the  dry  adiabatic  lapse  rate.  In  Figure  52  the  contributions  to  Cn2  from  the 
temperature  and  pressure  gradient  terms  are  plotted  individually,  along  with  the  resulting 
C2  including  both  terms. 


Figure  52.  Comparison  of  contributions  to  Cn  by  temperature  and  pressure 
gradient  terms.  The  right  plot  is  the  lowest  12  km  of  the  same  data  as  the  left  plot, 
expanded  to  show  detail.  This  Dayton,  OH  profile  is  for  10  Nov,  2010,  0000  UTC. 

In  Figure  52,  the  pressure  gradient  term  is  much  smaller  than  the  temperature 
gradient  term  (by  about  two  orders  of  magnitude)  throughout  most  of  the  profile,  but  at  the 
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heights  where  the  temperature  gradient  term  drops  sharply  due  to  the  near  constant 
potential  temperature,  the  pressure  gradient  term  becomes  significant.  To  analyze  the 
improvement  this  modification  has  on  derived  Cn2  values,  a  series  of  comparisons  are  made 
to  scintillometer-based  measurements  made  by  Dr.  Mikhail  Vorontsov  and  his  research 
team  at  the  University  of  Dayton  [66].  These  measurements  are  from  November  and 
December  of  2010,  along  a  7  km  path  between  the  University  of  Dayton  and  the  Veterans 
Affairs  Medical  Center,  shown  in  Figure  53. 


Figure  53.  Depiction  of  the  7  km  University  of  Dayton  scintillometer  path  and 
nearest  ECMWF  grid  point.  Map  data  and  imagery:  Google  Earth,  Landsat  [66] 
[67]. 

These  two  months  of  University  of  Dayton  scintillometer  measurements  are  during 
the  time  period  covered  by  the  ECMWF  ERA-20C  data.  For  these  Cn2  calculations,  the 
values  of  Kh/Km  =1  and  Lo=  10m  are  used  since  they  are  more  appropriate  values  for  the 
boundary  layer,  and  the  constant  Kh/Km  prevents  any  problems  being  introduced  by  the 
gradient  Richardson  number  calculation.  The  ECMWF  data  are  for  a  single  grid  point  near 
the  scintillometer  path,  just  north  of  the  path  toward  the  VA  Medical  Center  end  of  the 
path. 
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The  Cn  values  (temperature  gradient  only,  and  combined  temperature  and  pressure 
gradients)  calculated  for  the  lowest  ECMWF  level  at  the  grid  point  shown  in  Figure  53  are 
plotted  in  Figure  54,  along  with  the  UD  scintillometer  measured  Cn  values  along  the  7  km 
path.  These  data  represent  all  available  scintillometer  data  for  this  5-week  time  period, 
with  gaps  when  no  data  were  collected.  Figure  55  and  Figure  56  show  two  6-day  subsets 
of  these  data. 


Figure  54.  Analysis  of  modified  and  unmodified  Cn  derivation  technique,  each 
compared  to  UD  scintillometer  data  collected  by  the  University  of  Dayton  during 
Nov-Dec  2010.  For  all  three  of  these  plots,  the  dates/times  on  the  horizontal  axis  are 
UTC,  and  the  shaded  regions  identify  the  time  periods  between  local  sunset  and 
sunrise  [66]. 
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Figure  55.  Subset  of  Figure  54  data.  11  Nov  2010-17  Nov  2010  [66]. 


Figure  56.  Subset  of  Figure  54  data.  28  Nov  2010-4  Dec  2010  [66]. 

Throughout  this  range  of  dates,  it  is  observed  that  when  the  lower  atmosphere  is 
well-mixed  and  neutrally  buoyant,  the  temperature-only  C2  value  drops  by  one  to  three 
orders  of  magnitude.  At  these  low  values,  the  pressure  term  has  the  potential  to  become 
the  dominant  term,  leading  to  Cn2  values  closer  those  measured  by  the  scintillometer.  For 
the  33  days  that  scintillometer  data  were  available,  the  RMS  error  for  the  unmodified 
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technique  result  is  0.71  orders  of  magnitude.  With  the  pressure  gradient  term  added,  the 
RMS  error  is  reduced  to  0.65  orders  of  magnitude.  The  most  significant  improvement  is 
seen  during  the  quiescent  periods  in  turbulent  activity,  conditions  that  are  of  particular 
interest  to  directed  energy  research  efforts  as  their  favorable  conditions  are  often  targeted 
for  system  testing  and  evaluation. 
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5  Conclusions 


5.1  Summary  of  Findings 

This  research  develops  a  method  to  exploit  AIRS  temperature  retrievals  to  obtain  wind 
fields  and  other  atmospheric  parameters  ( Cf 2  and  C,,2)  of  importance  for  the  directed  energy 
community.  AIRS  data  takes  advantage  of  the  4.3  pm  carbon  dioxide  absorption  feature 
to  obtain  excellent  vertical  resolution  and  accurate  temperature  retrieval.  These  techniques 
are  not  documented  elsewhere,  and  can  be  applied  to  temperature  soundings  from  other 
satellite-based  infrared  sounders. 

One  of  the  conclusions  of  this  research  is  that,  even  though  the  AIRS  windfinding 
technique  has  some  proven  skill  to  determine  trends  in  wind  speed  and  direction,  it  does 
not  match  the  performance  of  modem  NWP  models.  The  volume  of  sounder  temperatures 
from  several  satellite-based  instruments  is  assimilated  into  the  weather  models,  so  the 
AIRS -based  technique  does  not  start  with  any  information  not  included  in  the  NWP 
models.  The  reason  the  NWP  winds  are  more  reliable  is  that  NWP  models  are  updated 
with  more  than  just  temperature  values.  Many  additional  sources  of  data,  such  as  surface 
and  upper-air  observations,  are  also  assimilated  by  the  model,  the  dynamics  of  which 
eliminate  many  of  the  measurement  errors  from  individual  observational  data  points.  After 
the  comparisons  detailed  in  this  document  it  is  clear  that,  even  when  remotely-sensed 
observations  of  the  needed  atmospheric  parameter  are  available,  NWP  data  should  be 
utilized. 

An  application  for  which  the  AIRS  windfinding  technique  is  still  well-suited  is 
generating  first  guess  winds  at  altitudes  where  NWP  model  data  is  not  available  or  easily 
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accessed  and  no  observations  are  available.  It  is  currently  difficult  (or  expensive)  to  access 
forecast  data  for  the  upper  stratosphere,  and  since  rocketsonde  campaigns  are  infrequent, 
direct  wind  observations  are  rare.  The  AIRS  windfinding  technique,  building  wind  profiles 
upward  from  the  top  pressure  level  of  GFS  or  another  NWP  model,  can  provide  some 
insight  into  wind  speed  and  direction  at  altitudes  where  no  other  information  is  easily 
accessible. 

The  modification  to  the  gradient  Richardson  number  is  not  found  to  be  of  significant 
benefit  to  estimating  Cn2.  In  atmospheric  layers  where  the  horizontal  wind  components  do 
not  change  between  the  bottom  of  the  layer  to  the  top,  the  inclusion  of  the  additional  wind 
shear  components  leads  to  a  decrease  in  Richardson  number  and  a  corresponding  increase 
in  calculated  Cn2.  The  additional  terms  are  very  small,  however,  and  the  modification  does 
not  lead  to  a  meaningful  increase  in  the  final  Cn2  value. 

A  significant  technique  developed  in  this  research  is  the  process  to  exploit  a  non¬ 
hydrostatic  model,  using  local  deviations  from  hydrostatic  equilibrium  to  estimate  the  size 
of  the  transient  pressure  fluctuations.  This  demonstrates  the  capability  to  calculate  and 
make  use  of  a  source  term  for  optical  turbulence  that  had  been  previously  neglected.  The 
small,  but  measurable  increase  in  Cn2  becomes  significant  in  a  well-mixed  layer,  where  the 
traditional  Tatarskii  approach  is  likely  to  underestimate  the  strength  of  the  optical 
turbulence. 

5.2  Impact  on  Directed  Energy  Research  Community 

The  atmospheric  characterization  techniques  and  the  refinement  of  existing  index  of 
refraction  structure  function  derivation  methods  can  benefit  the  directed  energy  research 
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community.  Enhanced  knowledge  about  the  strength  of  optical  turbulence  and  extinction 
along  a  path  of  interest  will  improve  our  understanding  of  how  a  directed  energy  system  or 
sensor  will  perform.  The  ability  to  better  characterize  the  atmosphere  will  allow  better, 
more  accurate  tests  of  systems  as  they  are  being  developed  and  evaluated.  Through  the 
combination  of  near  real-time  satellite  measurements  and  numerical  weather  model  output, 
the  research  paves  the  way  to  a  critical  predictive  capability,  arming  a  user  with  more 
knowledge  of  how  an  operational  system  would  perform  for  a  specified  path  and  time. 
These  techniques  can  allow  an  operator  a  passive  means  to  know  the  expected  effective 
range  of  a  laser  weapon  system,  as  determined  by  the  current  state  of  the  atmosphere. 
Optical  turbulence  and  extinction  vary  by  path,  so  the  ability  to  compare  system 
performance  for  several  engagement  geometries  would  provide  a  tactical  advantage  and 
enable  an  informed  decision  prior  to  system  employment. 

The  ability  to  passively  obtain  expected  bounds  of  C2  values  for  a  path  of  interest 
would  be  valuable  to  an  operational  high-altitude  HEL  system.  Actual  characterization  of 
the  atmosphere  at  the  time  of  the  engagement  would  be  performed  by  a  laser  beacon,  but 
prior  to  the  engagement,  the  AIRS/NWP-based  characterization  could  inform  decisions  on 
best  altitude  for  engagement  or  optimum  engagement  order  for  multiple  targets.  In  the 
possible  case  where  it  is  necessary  to  engage  a  target  that  is  beyond  the  range  of  the  beacon 
(but  not  beyond  the  range  of  the  HEL)  the  AIRS/NWP-derived  Cn2  values  would  provide 
characterization  of  the  strength  of  the  optical  turbulence  for  the  far  end  of  the  path, 
providing  an  estimate  of  required  dwell  time  for  the  engagement. 
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5.3  Recommendations  for  Future  Work 


The  techniques  developed  to  generate  standard  Cn 2  models  for  polar  and  tropical 
maritime  locations  could  be  used  to  create  additional,  more  specialized  standard  models. 
By  narrowing  the  geographic  region,  or  selecting  only  soundings  during  a  specific  time  of 
day  or  time  of  year,  the  envelope  of  expected  Cn2  values  for  a  more  narrow  set  of  conditions 
could  be  established.  The  use  of  AIRS  data  limits  time  of  day  investigation  to  the  overpass 
time  (around  0130  and  1330  local  time)  but  by  using  NWP  model  data,  standard  models 
could  be  created  for  any  time  of  day  (limited  only  by  the  temporal  resolution  of  the  model 
selected). 

Another  interesting  area  for  future  research  with  the  AIRS  data  is  an  extension  of 
the  temperature  accuracy  study  detailed  in  Section  4. 1 .  For  the  two  Alaska  locations,  the 
study  found  that  the  temperature  error  above  100  hPa  was  greater  for  cloud-free  soundings 
than  soundings  with  high  cloud  fraction.  This  may  indicate  a  problem  in  the  AIRS 
algorithm  detecting  high-altitude  clouds.  Additional  study,  including  more  locations,  and 
making  comparison  to  cloud  cover  observations  from  other  sources  could  determine  the 
source  of  this  unexpected  result. 

There  are  opportunities  to  continue  the  pressure  gradient  research,  coupling  it  with 
AFIT’s  on-going  radar  research.  The  ability  to  directly  measure  pressure  gradients  with 
the  radar  while  simultaneously  calculating  the  pressure  fluctuations  from  a  non-hydrostatic 

(  K  /  4/ 

weather  model  would  allow  calculation  of  the  group  of  coefficients  a 2  H/v  L0/3 

V  / 

which  prove  very  difficult  to  derive  from  easily  measured  parameters.  This  would  require 
access  to  current  non-hydrostatic  model  data,  and  a  good  choice  would  be  the  WRF  model. 
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Running  a  customized  WRF  model  would  allow  control  of  the  model  resolution  and  better 
study  of  the  effects  of  changes  to  initialization  data,  but  will  come  with  the  computational 
expense  of  running  a  high-resolution,  non-hydrostatic  NWP  model  and  the  corresponding 
lengthy  processing  times. 
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